Public Member Functions | Private Member Functions | Private Attributes | List of all members
dunefd::IniSegReco Class Reference
Inheritance diagram for dunefd::IniSegReco:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Member Functions

 IniSegReco (fhicl::ParameterSet const &p)
 
 IniSegReco (IniSegReco const &)=delete
 
 IniSegReco (IniSegReco &&)=delete
 
IniSegRecooperator= (IniSegReco const &)=delete
 
IniSegRecooperator= (IniSegReco &&)=delete
 
void beginJob () override
 
void reconfigure (fhicl::ParameterSet const &p)
 
void produce (art::Event &e) override
 
- 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
 
Produceroperator= (Producer const &)=delete
 
Produceroperator= (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
 
Modifieroperator= (Modifier const &)=delete
 
Modifieroperator= (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

void ResetVars (void)
 
void chargeParticlesatVtx (art::Event const &evt)
 
bool insideFidVol (TLorentzVector const &pvtx) const
 
float t0Corr (art::Event const &evt, detinfo::DetectorPropertiesData const &detProp, TLorentzVector const &pvtx)
 
TVector2 getMCvtx2d (TVector3 const &mcvtx3d, const size_t cryo, const size_t tpc, const size_t plane) const
 
TVector2 getMCdir2d (TVector3 const &mcvtx3d, TVector3 const &mcdir3d, const size_t cryo, const size_t tpc, const size_t plane) const
 
TVector3 findElDir (art::Ptr< simb::MCTruth > const mctruth) const
 
std::vector< TVector3 > findPhDir () const
 
std::vector< TVector3 > findDirs (art::Ptr< simb::MCTruth > const mctruth, int pdg) const
 
void collectCls (art::Event const &evt, detinfo::DetectorPropertiesData const &detProp, art::Ptr< simb::MCTruth > const mctruth)
 
std::vector< dunefd::Hit2DreselectCls (std::map< size_t, std::vector< dunefd::Hit2D > > const &cls, art::Ptr< simb::MCTruth > const mctruth, const size_t cryo, const size_t tpc, const size_t plane) const
 
void make3dseg (art::Event const &evt, detinfo::DetectorPropertiesData const &detProp, std::vector< std::vector< Hit2D > > const &src, TVector3 const &primary)
 
recob::Track convertFrom (pma::Track3D const &src)
 
void UseClusters (art::Event const &evt, detinfo::DetectorPropertiesData const &detProp)
 
void UseTracks (art::Event const &evt, detinfo::DetectorPropertiesData const &detProp)
 

Private Attributes

std::vector< pma::Track3D * > pmatracks
 
TTree * fTree
 
int run
 
int subrun
 
int event
 
short isdata
 
int trkindex
 
Float_t lep_dedx
 
Float_t dx
 
Float_t lep_dist
 
Float_t lep_distx
 
Float_t lep_disty
 
Float_t lep_distz
 
Float_t lep_distreco
 
Float_t lep_distrecox
 
Float_t lep_distrecoy
 
Float_t lep_distrecoz
 
Float_t cos
 
Float_t coslx
 
Float_t cosly
 
Float_t coslz
 
Float_t coslrx
 
Float_t coslry
 
Float_t coslrz
 
Float_t t0
 
Float_t vtxrecomc
 
Float_t vtxrecomcx
 
Float_t vtxrecomcy
 
Float_t vtxrecomcz
 
Float_t vtxupstream
 
Float_t vtxupstreamx
 
Float_t vtxupstreamy
 
Float_t vtxupstreamz
 
int ngamma
 
Float_t convdist
 
std::string fHitsModuleLabel
 
std::string fClusterModuleLabel
 
std::string fVertexModuleLabel
 
std::string fTrackModuleLabel
 
std::string fGenieGenModuleLabel
 
double fFidVolCut
 
pma::ProjectionMatchingAlg fProjectionMatchingAlg
 
std::ofstream file
 
std::ofstream file1
 

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
ConsumesCollectorconsumesCollector ()
 
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 50 of file IniSegReco_module.cc.

Constructor & Destructor Documentation

dunefd::IniSegReco::IniSegReco ( fhicl::ParameterSet const &  p)
explicit

Definition at line 151 of file IniSegReco_module.cc.

152 : EDProducer{pset}, fProjectionMatchingAlg(pset.get< fhicl::ParameterSet >("ProjectionMatchingAlg"))
153 {
154  this->reconfigure(pset);
155  produces< std::vector<recob::Track> >();
156  produces< std::vector<recob::SpacePoint> >();
157  produces< art::Assns<recob::Track, recob::Hit> >();
158  produces< art::Assns<recob::Track, recob::SpacePoint> >();
159  produces< art::Assns<recob::SpacePoint, recob::Hit> >();
160 }
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
void reconfigure(fhicl::ParameterSet const &p)
pma::ProjectionMatchingAlg fProjectionMatchingAlg
dunefd::IniSegReco::IniSegReco ( IniSegReco const &  )
delete
dunefd::IniSegReco::IniSegReco ( IniSegReco &&  )
delete

Member Function Documentation

void dunefd::IniSegReco::beginJob ( )
overridevirtual

Reimplemented from art::EDProducer.

Definition at line 163 of file IniSegReco_module.cc.

164 {
165  // Implementation of optional member function here.
167  fTree = tfs->make<TTree>("nueana","analysis tree");
168  fTree->Branch("run",&run,"run/I");
169  fTree->Branch("subrun",&subrun,"subrun/I");
170  fTree->Branch("event",&event,"event/I");
171  fTree->Branch("lep_dedx",&lep_dedx,"lep_dedx/F");
172  fTree->Branch("dx",&dx,"dx/F");
173  fTree->Branch("lep_dist",&lep_dist,"lep_dist/F");
174  fTree->Branch("lep_distx",&lep_distx,"lep_distx/F");
175  fTree->Branch("lep_disty",&lep_disty,"lep_disty/F");
176  fTree->Branch("lep_distz",&lep_distz,"lep_distz/F");
177  fTree->Branch("lep_distreco",&lep_distreco,"lep_distreco/F");
178  fTree->Branch("lep_distrecox",&lep_distrecox,"lep_distrecox/F");
179  fTree->Branch("lep_distrecoy",&lep_distrecoy,"lep_distrecoy/F");
180  fTree->Branch("lep_distrecoz",&lep_distrecoz,"lep_distrecoz/F");
181  fTree->Branch("cos", &cos, "cos/F");
182  fTree->Branch("coslx", &coslx, "coslx/F");
183  fTree->Branch("cosly", &cosly, "cosly/F");
184  fTree->Branch("coslz", &coslz, "coslz/F");
185  fTree->Branch("coslrx", &coslrx, "coslrx/F");
186  fTree->Branch("coslry", &coslry, "coslry/F");
187  fTree->Branch("coslrz", &coslrz, "coslrz/F");
188  fTree->Branch("t0", &t0, "t0/F");
189  fTree->Branch("vtxrecomc", &vtxrecomc, "vtxrecomc/F");
190  fTree->Branch("vtxrecomcx", &vtxrecomcx, "vtxrecomcx/F");
191  fTree->Branch("vtxrecomcy", &vtxrecomcy, "vtxrecomcy/F");
192  fTree->Branch("vtxrecomcz", &vtxrecomcz, "vtxrecomcz/F");
193  fTree->Branch("vtxupstream", &vtxupstream, "vtxupstream/F");
194  fTree->Branch("vtxupstreamx", &vtxupstreamx, "vtxupstreamx/F");
195  fTree->Branch("vtxupstreamy", &vtxupstreamy, "vtxupstreamy/F");
196  fTree->Branch("vtxupstreamz", &vtxupstreamz, "vtxupstreamz/F");
197  fTree->Branch("ngamma", &ngamma, "ngamma/I");
198  fTree->Branch("convdist", &convdist, "convdist/F");
199 
200 
201  file.open("data.dat");
202  file1.open("data1.dat");
203 }
code to link reconstructed objects back to the MC truth information
Event finding and building.
void dunefd::IniSegReco::chargeParticlesatVtx ( art::Event const &  evt)
private

Definition at line 819 of file IniSegReco_module.cc.

820 {
821  std::cout << " chargeParticleatVtx " << std::endl;
823  const sim::ParticleList& plist = pi_serv->ParticleList();
824 
825  TVector3 pospri; bool primary = false;
826  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar)
827  {
828  const simb::MCParticle* particle = ipar->second;
829 
830  if (particle->Process() != "primary") continue;
831 
832  pospri.SetXYZ(particle->Position(0).X(), particle->Position(0).Y(), particle->Position(0).Z());
833  primary = true; break;
834  }
835 
836  if (!primary) return;
837 
838  size_t ninteractions = 0;
839  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar)
840  {
841  std::vector< double > temp;
842  const simb::MCParticle* particle = ipar->second;
843 
844  std::cout << " pdg " << particle->PdgCode() << " " << particle->P(0) << std::endl;
845  std::cout << " position " << particle->Position(0).X() << ", " << particle->Position(0).Y() << ", " << particle->Position(0).Z() << std::endl;
846 
847  const TLorentzVector& startpos = particle->Position(0);
848  TVector3 start(startpos.X(), startpos.Y(), startpos.Z());
849  const TLorentzVector& endpos = particle->EndPosition();
850  TVector3 stop(endpos.X(), endpos.Y(), endpos.Z());
851 
852  double diststartpri = std::sqrt(pma::Dist2(start, pospri));
853  double diststoppri = std::sqrt(pma::Dist2(stop, pospri));
854  if ((diststartpri > 0.1) && (diststartpri < 1.0))
855  {
856  //if (diststoppri < 3.0)
857  //{
858  double ek = particle->E(0) - particle->Mass();
859  file << evt.run() << " " << evt.id().event() << " " << particle->PdgCode() << " " << diststartpri << " " << diststoppri << " " << particle->P(0) << " " << particle->Mass() << " " << ek << std::endl;
860 
861  if (ek > 0.05) ninteractions++;
862  //}
863  }
864  }
865 
866  file1 << evt.run() << " " << evt.id().event() << " " << ninteractions << std::endl;
867 }
double E(const int i=0) const
Definition: MCParticle.h:233
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:219
int PdgCode() const
Definition: MCParticle.h:212
unsigned int event
Definition: DataStructs.h:636
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:225
unsigned int run
Definition: DataStructs.h:637
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:37
double Mass() const
Definition: MCParticle.h:239
intermediate_table::const_iterator const_iterator
std::string Process() const
Definition: MCParticle.h:215
double P(const int i=0) const
Definition: MCParticle.h:234
const sim::ParticleList & ParticleList() const
TCEvent evt
Definition: DataStructs.cxx:7
QTextStream & endl(QTextStream &s)
void dunefd::IniSegReco::collectCls ( art::Event const &  evt,
detinfo::DetectorPropertiesData const &  detProp,
art::Ptr< simb::MCTruth > const  mctruth 
)
private

Definition at line 657 of file IniSegReco_module.cc.

660 {
662 
663  // * clusters
664  std::vector<art::Ptr<recob::Cluster> > clusterlist;
665  auto clusterListHandle = evt.getHandle< std::vector<recob::Cluster> >(fClusterModuleLabel);
666  if (clusterListHandle)
667  {
668  art::fill_ptr_vector(clusterlist, clusterListHandle);
669  art::FindManyP< recob::Hit > hc(clusterListHandle, evt, fClusterModuleLabel);
670 
671  for (size_t c = 0; c < geom->Ncryostats(); ++c)
672  {
673  const geo::CryostatGeo& cryo = geom->Cryostat(c);
674  for (size_t t = 0; t < cryo.NTPC(); ++t)
675  {
676  const geo::TPCGeo& tpc = cryo.TPC(t);
677  std::vector< std::vector< dunefd::Hit2D > > clinput;
678  for (size_t p = 0; p < tpc.Nplanes(); ++p)
679  {
680  std::map<size_t, std::vector< dunefd::Hit2D > > cls;
681  for (size_t i = 0; i < clusterlist.size(); ++i)
682  {
683  std::vector< art::Ptr<recob::Hit> > hitscl;
684  std::vector< Hit2D > hits2dcl;
685  hitscl = hc.at(i);
686 
687  bool found = false;
688  for (size_t h = 0; h < hitscl.size(); ++h)
689  {
690  size_t wire = hitscl[h]->WireID().Wire;
691  size_t plane = hitscl[h]->WireID().Plane;
692  size_t tpc = hitscl[h]->WireID().TPC;
693  size_t cryo = hitscl[h]->WireID().Cryostat;
694 
695  if ((plane == p) && (tpc == t) && (cryo == c))
696  {
697  found = true;
698  TVector2 point = pma::WireDriftToCm(detProp, wire, hitscl[h]->PeakTime(), plane, tpc, cryo);
699  Hit2D hit2d(point, hitscl[h].key());
700  hits2dcl.push_back(hit2d);
701  }
702  }
703 
704  if (found) cls[i] = hits2dcl;
705  }
706 
707  if (cls.size())
708  clinput.push_back(reselectCls(cls, mctruth, c, t, p));
709  }
710 
711  if (clinput.size() > 1)
712  {
713  const TLorentzVector& pvtx = mctruth->GetNeutrino().Nu().Position();
714  TVector3 primary(pvtx.X(), pvtx.Y(), pvtx.Z());
715  make3dseg(evt, detProp, clinput, primary);
716  }
717  }
718  }
719 
720  }
721 }
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:219
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
void make3dseg(art::Event const &evt, detinfo::DetectorPropertiesData const &detProp, std::vector< std::vector< Hit2D > > const &src, TVector3 const &primary)
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:165
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
Geometry information for a single TPC.
Definition: TPCGeo.h:38
TVector2 WireDriftToCm(detinfo::DetectorPropertiesData const &detProp, unsigned int wire, float drift, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:294
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
def key(type, name=None)
Definition: graph.py:13
p
Definition: test.py:223
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
std::string fClusterModuleLabel
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
std::vector< dunefd::Hit2D > reselectCls(std::map< size_t, std::vector< dunefd::Hit2D > > const &cls, art::Ptr< simb::MCTruth > const mctruth, const size_t cryo, const size_t tpc, const size_t plane) const
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
recob::Track dunefd::IniSegReco::convertFrom ( pma::Track3D const &  src)
private

Definition at line 1046 of file IniSegReco_module.cc.

1047 {
1048  std::vector< TVector3 > xyz, dircos;
1049 
1050  for (size_t i = 0; i < src.size(); ++i)
1051  {
1052  xyz.push_back(src[i]->Point3D());
1053 
1054  if (i < src.size() - 1)
1055  {
1056  TVector3 dc(src[i + 1]->Point3D());
1057  dc -= src[i]->Point3D();
1058  dc *= 1.0 / dc.Mag();
1059  dircos.push_back(dc);
1060  }
1061  else dircos.push_back(dircos.back());
1062  }
1063 
1064  if (xyz.size() != dircos.size())
1065  {
1066  mf::LogError("IniSegReco") << "pma::Track3D to recob::Track conversion problem.";
1067  }
1070  recob::Track::Flags_t(xyz.size()), false),
1072 }
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 > > SMatrixSym55
Definition: TrackingTypes.h:85
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Definition: TrackingTypes.h:58
::fhicl::TupleAs< Point(::geo::Length_t,::geo::Length_t,::geo::Length_t)> Point3D
Atom object for reading a 3D point or vector (centimeters).
TrackTrajectory::Flags_t Flags_t
Definition: Track.h:68
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
A trajectory in space reconstructed from hits.
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:55
TrackCollectionProxyElement< TrackCollProxy > Track
Proxy to an element of a proxy collection of recob::Track objects.
Definition: Track.h:1036
std::vector< TVector3 > dunefd::IniSegReco::findDirs ( art::Ptr< simb::MCTruth > const  mctruth,
int  pdg 
) const
private

Definition at line 923 of file IniSegReco_module.cc.

924 {
925  std::vector< TVector3 > dirs;
926 
928  const sim::ParticleList& plist = pi_serv->ParticleList();
929 
930  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar)
931  {
932  simb::MCParticle* particle = ipar->second;
933  TLorentzVector mom = particle->Momentum();
934  TVector3 momvec(mom.Px(), mom.Py(), mom.Pz());
935 
936  if ((particle->PdgCode() == pdg) && (momvec.Mag() > 0.030))
937  {
938  TLorentzVector momconv = particle->EndMomentum();
939  TVector3 momconvec3(momconv.Px(), momconv.Py(), momconv.Pz());
940  TVector3 dir = momconvec3 * (1 / momconvec3.Mag());
941  dirs.push_back(dir);
942  }
943  }
944 
945  return dirs;
946 }
int PdgCode() const
Definition: MCParticle.h:212
intermediate_table::const_iterator const_iterator
string dir
const sim::ParticleList & ParticleList() const
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
const TLorentzVector & EndMomentum() const
Definition: MCParticle.h:240
TVector3 dunefd::IniSegReco::findElDir ( art::Ptr< simb::MCTruth > const  mctruth) const
private

Definition at line 802 of file IniSegReco_module.cc.

803 {
804  TVector3 dir(0, 0, 0);
805  const simb::MCParticle& lepton = mctruth->GetNeutrino().Lepton();
806 
807  if (abs(lepton.PdgCode()) == 11)
808  {
809  TLorentzVector mom = lepton.Momentum();
810  TVector3 momvec3(mom.Px(), mom.Py(), mom.Pz());
811  dir = momvec3 * (1 / momvec3.Mag());
812  }
813 
814  return dir;
815 }
int PdgCode() const
Definition: MCParticle.h:212
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
string dir
T abs(T value)
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
std::vector< TVector3 > dunefd::IniSegReco::findPhDir ( ) const
private

Definition at line 871 of file IniSegReco_module.cc.

872 {
873  std::vector< TVector3 > phdirs;
875  const sim::ParticleList& plist = pi_serv->ParticleList();
876 
877  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar)
878  {
879  const simb::MCParticle* particle = ipar->second;
880  if (particle->Process() != "primary") continue;
881 
882  TVector3 dir(0, 0, 0);
883  if (particle->PdgCode() == 111)
884  {
885  if (particle->NumberDaughters() != 2) continue;
886 
887  const simb::MCParticle* daughter1 = pi_serv->TrackIdToParticle_P(particle->Daughter(0));
888  if (daughter1->PdgCode() != 22) continue;
889 
890  const simb::MCParticle* daughter2 = pi_serv->TrackIdToParticle_P(particle->Daughter(1));
891  if (daughter2->PdgCode() != 22) continue;
892 
893  if (daughter1->EndProcess() == "conv")
894  {
895  TLorentzVector mom = pi_serv->TrackIdToParticle_P(particle->Daughter(0))->Momentum();
896  TVector3 momvec3(mom.Px(), mom.Py(), mom.Pz());
897  dir = momvec3 * (1 / momvec3.Mag());
898  phdirs.push_back(dir);
899  }
900 
901  if (daughter2->EndProcess() == "conv")
902  {
903  TLorentzVector mom = pi_serv->TrackIdToParticle_P(particle->Daughter(1))->Momentum();
904  TVector3 momvec3(mom.Px(), mom.Py(), mom.Pz());
905  dir = momvec3 * (1 / momvec3.Mag());
906  phdirs.push_back(dir);
907  }
908  }
909  else if (particle->PdgCode() == 22)
910  {
911  TLorentzVector mom = particle->Momentum();
912  TVector3 momvec3(mom.Px(), mom.Py(), mom.Pz());
913  dir = momvec3 * (1 / momvec3.Mag());
914  phdirs.push_back(dir);
915  }
916  }
917 
918  return phdirs;
919 }
int PdgCode() const
Definition: MCParticle.h:212
const simb::MCParticle * TrackIdToParticle_P(int id) const
intermediate_table::const_iterator const_iterator
string dir
std::string Process() const
Definition: MCParticle.h:215
int NumberDaughters() const
Definition: MCParticle.h:217
int Daughter(const int i) const
Definition: MCParticle.cxx:112
std::string EndProcess() const
Definition: MCParticle.h:216
const sim::ParticleList & ParticleList() const
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
TVector2 dunefd::IniSegReco::getMCdir2d ( TVector3 const &  mcvtx3d,
TVector3 const &  mcdir3d,
const size_t  cryo,
const size_t  tpc,
const size_t  plane 
) const
private

Definition at line 781 of file IniSegReco_module.cc.

783 {
784  TVector3 shift3d = mcvtx3d + mcdir3d;
785  TVector2 shift2d = pma::GetProjectionToPlane(shift3d, plane, tpc, cryo) - getMCvtx2d(mcvtx3d, cryo, tpc, plane);
786  TVector2 mcdir2d = shift2d * (1.0 / shift2d.Mod());
787 
788  return mcdir2d;
789 }
TVector2 getMCvtx2d(TVector3 const &mcvtx3d, const size_t cryo, const size_t tpc, const size_t plane) const
TVector2 GetProjectionToPlane(const TVector3 &p, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:270
TVector2 dunefd::IniSegReco::getMCvtx2d ( TVector3 const &  mcvtx3d,
const size_t  cryo,
const size_t  tpc,
const size_t  plane 
) const
private

Definition at line 793 of file IniSegReco_module.cc.

794 {
795  TVector2 mcvtx2d = pma::GetProjectionToPlane(mcvtx3d, plane, tpc, cryo);
796 
797  return mcvtx2d;
798 }
TVector2 GetProjectionToPlane(const TVector3 &p, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:270
bool dunefd::IniSegReco::insideFidVol ( TLorentzVector const &  pvtx) const
private

Definition at line 986 of file IniSegReco_module.cc.

987 {
989  double vtx[3] = {pvtx.X(), pvtx.Y(), pvtx.Z()};
990  bool inside = false;
991 
992  geo::TPCID idtpc = geom->FindTPCAtPosition(vtx);
993 
994  if (geom->HasTPC(idtpc))
995  {
996  const geo::TPCGeo& tpcgeo = geom->GetElement(idtpc);
997  double minx = tpcgeo.MinX(); double maxx = tpcgeo.MaxX();
998  double miny = tpcgeo.MinY(); double maxy = tpcgeo.MaxY();
999  double minz = tpcgeo.MinZ(); double maxz = tpcgeo.MaxZ();
1000 
1001  /*for (size_t c = 0; c < geom->Ncryostats(); c++)
1002  {
1003  const geo::CryostatGeo& cryostat = geom->Cryostat(c);
1004  for (size_t t = 0; t < cryostat.NTPC(); t++)
1005  {
1006  const geo::TPCGeo& tpcg = cryostat.TPC(t);
1007  if (tpcg.MinX() < minx) minx = tpcg.MinX();
1008  if (tpcg.MaxX() > maxx) maxx = tpcg.MaxX();
1009  if (tpcg.MinY() < miny) miny = tpcg.MinY();
1010  if (tpcg.MaxY() > maxy) maxy = tpcg.MaxY();
1011  if (tpcg.MinZ() < minz) minz = tpcg.MinZ();
1012  if (tpcg.MaxZ() > maxz) maxz = tpcg.MaxZ();
1013  }
1014  }*/
1015 
1016  //x
1017  double dista = fabs(minx - pvtx.X());
1018  double distb = fabs(pvtx.X() - maxx);
1019 ;
1020  if ((pvtx.X() > minx) && (pvtx.X() < maxx) &&
1021  (dista > fFidVolCut) && (distb > fFidVolCut))
1022  {
1023  inside = true;
1024  }
1025  else { inside = false; }
1026  //y
1027  dista = fabs(maxy - pvtx.Y());
1028  distb = fabs(pvtx.Y() - miny);
1029  if (inside && (pvtx.Y() > miny) && (pvtx.Y() < maxy) &&
1030  (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
1031  else inside = false;
1032 
1033  //z
1034  dista = fabs(maxz - pvtx.Z());
1035  distb = fabs(pvtx.Z() - minz);
1036  if (inside && (pvtx.Z() > minz) && (pvtx.Z() < maxz) &&
1037  (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
1038  else inside = false;
1039  }
1040 
1041  return inside;
1042 }
CryostatGeo const & GetElement(geo::CryostatID const &cryoid) const
double MinX() const
Returns the world x coordinate of the start of the box.
Definition: BoxBoundedGeo.h:88
Geometry information for a single TPC.
Definition: TPCGeo.h:38
double MaxX() const
Returns the world x coordinate of the end of the box.
Definition: BoxBoundedGeo.h:91
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
double MinZ() const
Returns the world z coordinate of the start of the box.
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
double MaxY() const
Returns the world y coordinate of the end of the box.
bool HasTPC(geo::TPCID const &tpcid) const
Returns whether we have the specified TPC.
double MaxZ() const
Returns the world z coordinate of the end of the box.
double MinY() const
Returns the world y coordinate of the start of the box.
void dunefd::IniSegReco::make3dseg ( art::Event const &  evt,
detinfo::DetectorPropertiesData const &  detProp,
std::vector< std::vector< Hit2D > > const &  src,
TVector3 const &  primary 
)
private

Definition at line 752 of file IniSegReco_module.cc.

757 {
758  if (src.size() < 2) return;
759  if ((src[0].size() < 2) || (src[1].size() < 2)) return;
760 
761  std::vector<art::Ptr<recob::Hit> > hitlist;
762  auto hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
763  if (hitListHandle)
764  art::fill_ptr_vector(hitlist, hitListHandle);
765 
766  std::vector< art::Ptr<recob::Hit> > hitsel;
767 
768  for (size_t i = 0; i < src.size(); ++i)
769  for (size_t h = 0; h < src[i].size(); ++h)
770  hitsel.push_back(hitlist[src[i][h].GetKey()]);
771 
772  if (hitsel.size() > 5)
773  {
775  pmatracks.push_back(trk);
776  }
777 }
std::vector< pma::Track3D * > pmatracks
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
pma::Track3D * buildSegment(const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits_1, const std::vector< art::Ptr< recob::Hit >> &hits_2={}) const
std::string fHitsModuleLabel
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
pma::ProjectionMatchingAlg fProjectionMatchingAlg
IniSegReco& dunefd::IniSegReco::operator= ( IniSegReco const &  )
delete
IniSegReco& dunefd::IniSegReco::operator= ( IniSegReco &&  )
delete
void dunefd::IniSegReco::produce ( art::Event e)
overridevirtual

Implements art::EDProducer.

Definition at line 247 of file IniSegReco_module.cc.

248 {
249  ResetVars();
251 
252  run = evt.run();
253  subrun = evt.subRun();
254  event = evt.id().event();
255  isdata = evt.isRealData();
256 
257  std::unique_ptr< std::vector< recob::Track > > tracks(new std::vector< recob::Track >);
258  std::unique_ptr< std::vector< recob::SpacePoint > > allsp(new std::vector< recob::SpacePoint >);
259 
260  std::unique_ptr< art::Assns< recob::Track, recob::Hit > > trk2hit(new art::Assns< recob::Track, recob::Hit >);
261  std::unique_ptr< art::Assns< recob::Track, recob::SpacePoint > > trk2sp(new art::Assns< recob::Track, recob::SpacePoint >);
262  std::unique_ptr< art::Assns< recob::SpacePoint, recob::Hit > > sp2hit(new art::Assns< recob::SpacePoint, recob::Hit >);
263 
264  TVector3 primary(0, 0, 0);
265  bool isinside = false;
266  if (!isdata)
267  {
268  // * MC truth information
269  std::vector<art::Ptr<simb::MCTruth> > mclist;
270  auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
271  if (mctruthListHandle)
272  art::fill_ptr_vector(mclist, mctruthListHandle);
273 
274  if (mclist.size())
275  {
276  art::Ptr<simb::MCTruth> mctruth = mclist[0];
277  const TLorentzVector& pvtx = mctruth->GetNeutrino().Nu().Position();
278  primary = TVector3(pvtx.X(), pvtx.Y(), pvtx.Z());
279 
280  if (insideFidVol(pvtx))
281  {
282  isinside = true;
284 
286  const sim::ParticleList& plist = pi_serv->ParticleList();
287 
288  bool photon = false;
289  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar)
290  {
291  simb::MCParticle* particle = ipar->second;
292  TLorentzVector mom = particle->Momentum();
293  TVector3 momvec(mom.Px(), mom.Py(), mom.Pz());
294 
295  if ((particle->PdgCode() == 22) && (momvec.Mag() > 0.030))
296  {
297  ngamma++; photon = true;
298  TLorentzVector conversion = particle->EndPosition();
299  TVector3 convec(conversion.X(), conversion.Y(), conversion.Z());
300  convdist = std::sqrt(pma::Dist2(primary, convec));
301  }
302  }
303  if (!photon) ngamma = -10;
304  }
305  }
306  }
307 
308  if (isinside)
309  {
310  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt);
311  UseTracks(evt, detProp);
312  fTree->Fill();
313  }
314 
315  evt.put(std::move(tracks));
316  evt.put(std::move(allsp));
317  evt.put(std::move(trk2hit));
318  evt.put(std::move(trk2sp));
319  evt.put(std::move(sp2hit));
320 }
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:219
int PdgCode() const
Definition: MCParticle.h:212
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
unsigned int event
Definition: DataStructs.h:636
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:225
unsigned int run
Definition: DataStructs.h:637
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:37
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
void chargeParticlesatVtx(art::Event const &evt)
intermediate_table::const_iterator const_iterator
std::string fGenieGenModuleLabel
def move(depos, offset)
Definition: depos.py:107
bool insideFidVol(TLorentzVector const &pvtx) const
const sim::ParticleList & ParticleList() const
Definition: tracks.py:1
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
void UseTracks(art::Event const &evt, detinfo::DetectorPropertiesData const &detProp)
unsigned int subRun
Definition: DataStructs.h:638
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
void dunefd::IniSegReco::reconfigure ( fhicl::ParameterSet const &  p)

Definition at line 205 of file IniSegReco_module.cc.

206 {
207  fHitsModuleLabel = pset.get< std::string >("HitsModuleLabel");
208  fClusterModuleLabel = pset.get< std::string >("ClusterModuleLabel");
209  fTrackModuleLabel = pset.get< std::string >("TrackModuleLabel");
210  fVertexModuleLabel = pset.get< std::string >("VertexModuleLabel");
211  fGenieGenModuleLabel = pset.get< std::string >("GenieGenModuleLabel");
212 
213  fFidVolCut = pset.get< double >("FidVolCut");
214  return;
215 }
std::string string
Definition: nybbler.cc:12
std::string fGenieGenModuleLabel
std::string fHitsModuleLabel
std::string fVertexModuleLabel
std::string fClusterModuleLabel
std::string fTrackModuleLabel
std::vector< dunefd::Hit2D > dunefd::IniSegReco::reselectCls ( std::map< size_t, std::vector< dunefd::Hit2D > > const &  cls,
art::Ptr< simb::MCTruth > const  mctruth,
const size_t  cryo,
const size_t  tpc,
const size_t  plane 
) const
private

Definition at line 725 of file IniSegReco_module.cc.

726 {
727  std::vector< dunefd::Hit2D > cluster;
728  if (!cls.size()) return cluster;
729 
730  const simb::MCParticle& particle = mctruth->GetNeutrino().Nu();
731  const TLorentzVector& pvtx = particle.Position();
732 
733  if (insideFidVol(pvtx) && (abs(mctruth->GetNeutrino().Lepton().PdgCode()) == 11))
734  {
735  // mc
736  TVector3 mcvtx3d(pvtx.X(), pvtx.Y(), pvtx.Z());
737  TVector3 mcdir3d = findElDir(mctruth);
738  TVector2 mcvtx2d = getMCvtx2d(mcvtx3d, cryo, tpc, plane);
739  TVector2 mcdir2d = getMCdir2d(mcvtx3d, mcdir3d, cryo, tpc, plane);
740 
741  // reco: find the best clusters to proceed with segment reconstruction
742  IniSegAlg recoini(cls);
743  recoini.FeedwithMc(mcvtx2d, mcdir2d, mcdir3d);
744  cluster = recoini.GetCl();
745  }
746 
747  return cluster;
748 }
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:219
int PdgCode() const
Definition: MCParticle.h:212
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
void cluster(In first, In last, Out result, Pred *pred)
Definition: NNClusters.h:41
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
TVector2 getMCvtx2d(TVector3 const &mcvtx3d, const size_t cryo, const size_t tpc, const size_t plane) const
T abs(T value)
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
bool insideFidVol(TLorentzVector const &pvtx) const
TVector3 findElDir(art::Ptr< simb::MCTruth > const mctruth) const
TVector2 getMCdir2d(TVector3 const &mcvtx3d, TVector3 const &mcdir3d, const size_t cryo, const size_t tpc, const size_t plane) const
void dunefd::IniSegReco::ResetVars ( void  )
private

Definition at line 217 of file IniSegReco_module.cc.

218 {
219  ngamma = 0;
220  convdist = -10.0F;
221  pmatracks.clear();
222  lep_dist = -9999; //cm
223  lep_distx = -9999;
224  lep_disty = -9999;
225  lep_distz = -9999;
226  lep_distreco = -9999;
227  lep_distrecox = -9999;
228  lep_distrecoy = -9999;
229  lep_distrecoz = -9999;
230  lep_dedx = -9999;
231  t0 = -9999;
232  dx = -9999;
233  cos = -9999;
234  coslx = -9999; cosly = -9999; coslz = -9999;
235  coslrx = -9999; coslry = -9999; coslrz = -9999;
236  vtxrecomc = -9999;
237  vtxrecomcx = -9999;
238  vtxrecomcy = -9999;
239  vtxrecomcz = -9999;
240  vtxupstream = -9999;
241  vtxupstreamx = -9999;
242  vtxupstreamy = -9999;
243  vtxupstreamz = -9999;
244  return;
245 }
code to link reconstructed objects back to the MC truth information
std::vector< pma::Track3D * > pmatracks
float dunefd::IniSegReco::t0Corr ( art::Event const &  evt,
detinfo::DetectorPropertiesData const &  detProp,
TLorentzVector const &  pvtx 
)
private

Definition at line 950 of file IniSegReco_module.cc.

953 {
954  float corrt0x = 0.0F;
955 
957 
958  // * MC truth information
959  std::vector<art::Ptr<simb::MCTruth> > mclist;
960  auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
961  if (mctruthListHandle)
962  art::fill_ptr_vector(mclist, mctruthListHandle);
963 
964  double vtx[3] = {pvtx.X(), pvtx.Y(), pvtx.Z()};
965 
966  geo::TPCID idtpc = geom->FindTPCAtPosition(vtx);
967 
968  if (geom->HasTPC(idtpc))
969  if (mclist.size())
970  {
971  art::Ptr<simb::MCTruth> mctruth = mclist[0];
972  if (mctruth->NParticles())
973  {
974  simb::MCParticle particle = mctruth->GetParticle(0);
975  t0 = particle.T(); // ns
976  corrt0x = t0 * 1.e-3 * detProp.DriftVelocity();
977  }
978  }
979 
980 
981  return corrt0x;
982 }
code to link reconstructed objects back to the MC truth information
int NParticles() const
Definition: MCTruth.h:75
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
std::string fGenieGenModuleLabel
double T(const int i=0) const
Definition: MCParticle.h:224
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
bool HasTPC(geo::TPCID const &tpcid) const
Returns whether we have the specified TPC.
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
void dunefd::IniSegReco::UseClusters ( art::Event const &  evt,
detinfo::DetectorPropertiesData const &  detProp 
)
private

Definition at line 324 of file IniSegReco_module.cc.

326 {
327  if (!isdata)
328  {
329  // * MC truth information
330  std::vector<art::Ptr<simb::MCTruth> > mclist;
331  auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
332  if (mctruthListHandle)
333  art::fill_ptr_vector(mclist, mctruthListHandle);
334 
335  if (mclist.size())
336  {
337  art::Ptr<simb::MCTruth> mctruth = mclist[0];
338  collectCls(evt, detProp, mctruth);
339  }
340  }
341 }
std::string fGenieGenModuleLabel
void collectCls(art::Event const &evt, detinfo::DetectorPropertiesData const &detProp, art::Ptr< simb::MCTruth > const mctruth)
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
void dunefd::IniSegReco::UseTracks ( art::Event const &  evt,
detinfo::DetectorPropertiesData const &  detProp 
)
private

Definition at line 345 of file IniSegReco_module.cc.

347 {
348  if (!isdata)
349  {
350  // * hits
351  std::vector<art::Ptr<recob::Hit> > hitlist;
352  auto hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
353  if (hitListHandle)
354  art::fill_ptr_vector(hitlist, hitListHandle);
355 
356  // * tracks
357  std::vector<art::Ptr<recob::Track> > tracklist;
358  auto trackListHandle = evt.getHandle< std::vector<recob::Track> >(fTrackModuleLabel);
359  if (trackListHandle)
360  art::fill_ptr_vector(tracklist, trackListHandle);
361 
362  // * vertices
363  std::vector<art::Ptr<recob::Vertex> > vtxlist;
364  auto vtxListHandle = evt.getHandle< std::vector<recob::Vertex> >(fVertexModuleLabel);
365  if (vtxListHandle)
366  art::fill_ptr_vector(vtxlist, vtxListHandle);
367 
368  // * monte carlo
369  std::vector<art::Ptr<simb::MCTruth> > mclist;
370  auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
371  if (mctruthListHandle)
372  art::fill_ptr_vector(mclist, mctruthListHandle);
373 
374  if (mclist.size())
375  {
376  art::Ptr<simb::MCTruth> mctruth = mclist[0];
377 
378  const simb::MCParticle& particle = mctruth->GetNeutrino().Nu();
379  const TLorentzVector& pvtx = particle.Position();
380  TVector3 primary(pvtx.X(), pvtx.Y(), pvtx.Z());
381 
382  // search for the closest reco vertex to mc primary
383  TVector3 closestvtx; TVector3 minzvtx;
384  if (vtxlist.size())
385  {
386  double xyz[3] = {0.0, 0.0, 0.0};
387  vtxlist[0]->XYZ(xyz);
388  float vxreco = xyz[0];
389  if (vxreco > 0) vxreco -= t0Corr(evt, detProp, pvtx);
390  else vxreco += t0Corr(evt, detProp, pvtx);
391  xyz[0] = vxreco;
392 
393  TVector3 vtxreco(xyz);
394  vtxreco.SetXYZ(xyz[0], xyz[1], xyz[2]);
395  closestvtx.SetXYZ(xyz[0], xyz[1], xyz[2]);
396 
397  double mindist2 = pma::Dist2(primary, vtxreco);
398  // loop over vertices to look for the closest to primary
399  for (size_t v = 1; v < vtxlist.size(); ++v)
400  {
401  vtxlist[v]->XYZ(xyz);
402  float temp = xyz[0];
403  if (temp > 0) temp -= t0Corr(evt, detProp, pvtx);
404  else temp += t0Corr(evt, detProp, pvtx);
405  xyz[0] = temp;
406 
407  vtxreco.SetXYZ(xyz[0], xyz[1], xyz[2]);
408  float dist2 = pma::Dist2(primary, vtxreco);
409  if (dist2 < mindist2)
410  {
411  closestvtx.SetXYZ(xyz[0], xyz[1], xyz[2]);
412  mindist2 = dist2;
413  }
414  }
415 
416  // loop over vertices to look the most upstream
417  double minz_xyz[3] = {0.0, 0.0, 0.0};
418  double minz = 9999;
419  for (size_t v = 0; v < vtxlist.size(); ++v)
420  {
421  vtxlist[v]->XYZ(minz_xyz);
422  float temp = minz_xyz[0];
423  if (temp > 0) temp -= t0Corr(evt, detProp, pvtx);
424  else temp += t0Corr(evt, detProp, pvtx);
425  minz_xyz[0] = temp;
426 
427  if (minz_xyz[2] < minz)
428  {
429  minz = minz_xyz[2];
430  minzvtx.SetXYZ(minz_xyz[0], minz_xyz[1], minz_xyz[2]);
431  }
432  }
433 
434  // loop over tracks
435  for (size_t t = 0; t < tracklist.size(); ++t)
436  {
437  if (!tracklist[t]->NumberTrajectoryPoints()) continue;
438  TVector3 pos_p = tracklist[t]->LocationAtPoint<TVector3>(0);
439 
440  float temp = pos_p.X();
441  if (temp > 0) temp -= t0Corr(evt, detProp, pvtx);
442  else temp += t0Corr(evt, detProp, pvtx);
443  pos_p.SetX(temp);
444 
445  float dist2 = pma::Dist2(primary, pos_p);
446  if (dist2 < mindist2)
447  {
448  closestvtx.SetXYZ(pos_p.X(), pos_p.Y(), pos_p.Z());
449  mindist2 = dist2;
450  }
451 
452  if (pos_p.Z() < minzvtx.Z())
453  {
454  minzvtx.SetXYZ(pos_p.X(), pos_p.Y(), pos_p.Z());
455  }
456 
457  TVector3 pos_end = tracklist[t]->LocationAtPoint<TVector3>(tracklist[t]->NumberTrajectoryPoints()-1);
458  temp = pos_end.X();
459  if (temp > 0) temp -= t0Corr(evt, detProp, pvtx);
460  else temp += t0Corr(evt, detProp, pvtx);
461  pos_end.SetX(temp);
462 
463  dist2 = pma::Dist2(primary, pos_end);
464  if (dist2 < mindist2)
465  {
466  closestvtx.SetXYZ(pos_end.X(), pos_end.Y(), pos_end.Z());
467  mindist2 = dist2;
468  }
469 
470  if (pos_end.Z() < minzvtx.Z())
471  {
472  minzvtx.SetXYZ(pos_end.X(), pos_end.Y(), pos_end.Z());
473  }
474  }
475 
476  vtxrecomc = std::sqrt(pma::Dist2(primary, closestvtx));
477  vtxrecomcx = primary.X() - closestvtx.X();
478  vtxrecomcy = primary.Y() - closestvtx.Y();
479  vtxrecomcz = primary.Z() - closestvtx.Z();
480  vtxupstream = std::sqrt(pma::Dist2(primary, minzvtx));
481  vtxupstreamx = primary.X() - minzvtx.X();
482  vtxupstreamy = primary.Y() - minzvtx.Y();
483  vtxupstreamz = primary.Z() - minzvtx.Z();
484  }
485 
486  // vertex region of event
487  // background events
488  if (insideFidVol(pvtx) && (abs(mctruth->GetNeutrino().Lepton().PdgCode()) != 11) && tracklist.size())
489  {
490  // mc
491  TVector3 mcvtx3d(pvtx.X(), pvtx.Y(), pvtx.Z());
492  std::vector< TVector3 > mcdir3d = findPhDir();
493 
494  float minlep_dist = 9999;
495  for (size_t v = 0; v < mcdir3d.size(); ++v)
496  {
497  // reco
498  IniSegAlg recoini(tracklist, mcvtx3d);
499  recoini.FeedwithMc(mcdir3d[v]);
500 
501  if (recoini.IsFound())
502  {
503  art::Ptr<recob::Track> recotrack = recoini.GetTrk();
504  if (!recotrack->NumberTrajectoryPoints()) continue;
505 
506  art::FindManyP< recob::Hit > fb(trackListHandle, evt, fTrackModuleLabel);
507  std::vector< art::Ptr<recob::Hit> > recoinihit = fb.at(recotrack.key());
508 
509  // use recob::Track functionality as much as possible
510  // const double setlength = 2.5; double length = 0.0; // cm
511 
512  TVector3 pos_p = recotrack->LocationAtPoint<TVector3>(0);
513  float px = pos_p.X();
514  if (px > 0) px -= t0Corr(evt, detProp, pvtx);
515  else px += t0Corr(evt, detProp, pvtx);
516  pos_p.SetX(px);
517 
518  float ldist = std::sqrt(pma::Dist2(primary, pos_p));
519  if (ldist < minlep_dist)
520  {
521  minlep_dist = ldist;
522  lep_dist = ldist;
523  lep_distx = primary.X() - pos_p.X();
524  lep_disty = primary.Y() - pos_p.Y();
525  lep_distz = primary.Z() - pos_p.Z();
526  if (mclist.size())
527  {
528  lep_distreco = std::sqrt(pma::Dist2(closestvtx, pos_p));
529  lep_distrecox = closestvtx.X() - pos_p.X();
530  lep_distrecoy = closestvtx.Y() - pos_p.Y();
531  lep_distrecoz = closestvtx.Z() - pos_p.Z();
532  }
533 
534  cos = recoini.GetCos();
535  coslrx = mcdir3d[v].X(); coslry = mcdir3d[v].Y(); coslrz = mcdir3d[v].Z();
536 
537  /*************************************************************/
538  /* WARNING */
539  /*************************************************************/
540  /* The dQdx information in recob::Track has been deprecated */
541  /* since 2016 and in 11/2018 the recob::Track interface was */
542  /* changed and DQdxAtPoint and NumberdQdx were removed. */
543  /* Therefore the code below is now commented out */
544  /* (note that it was most likely not functional anyways). */
545  /* For any issue please contact: larsoft-team@fnal.gov */
546  /*************************************************************/
547  /*
548  size_t fp = 0; bool hitcoll = false;
549  for (size_t p = 0; p < recotrack->NumberTrajectoryPoints(); ++p)
550  if (recotrack->DQdxAtPoint(p, geo::kZ) > 0)
551  {pos_p = recotrack->LocationAtPoint<TVector3>(p); fp = p; hitcoll = true; break;}
552 
553  // loop over trajectory point to get dQdx.
554  if (hitcoll)
555  for (size_t p = (fp+1); p < recotrack->NumberTrajectoryPoints(); ++p)
556  {
557  TVector3 pos = recotrack->LocationAtPoint<TVector3>(p);
558  length += std::sqrt(pma::Dist2(pos_p, pos));
559  pos_p = recotrack->LocationAtPoint<TVector3>(p);
560 
561  if (length > setlength) break;
562  dx = length;
563  double dqdx_p = recotrack->DQdxAtPoint(p, geo::kZ);
564  if (dqdx_p > 0) lep_dedx += recoinihit[p]->SummedADC();
565  }
566 
567  if (dx > 0.0) lep_dedx /= dx;
568  */
569  /*************************************************************/
570  }
571  }
572  }
573  }
574  else if (insideFidVol(pvtx) && (abs(mctruth->GetNeutrino().Lepton().PdgCode()) == 11) && tracklist.size()) // from nue vertex
575  {
576  // mc
577  TVector3 mcvtx3d(pvtx.X(), pvtx.Y(), pvtx.Z());
578  TVector3 mcdir3d = findElDir(mctruth);
579  coslx = mcdir3d.X(); cosly = mcdir3d.Y(); coslz = mcdir3d.Z();
580 
581  // reco
582  IniSegAlg recoini(tracklist, mcvtx3d);
583  recoini.FeedwithMc(mcdir3d); // mcvtx3d == primary
584 
585  if (recoini.IsFound())
586  {
587  art::Ptr<recob::Track> recotrack = recoini.GetTrk();
588  if (!recotrack->NumberTrajectoryPoints()) return;
589 
590  art::FindManyP< recob::Hit > fb(trackListHandle, evt, fTrackModuleLabel);
591  std::vector< art::Ptr<recob::Hit> > recoinihit = fb.at(recotrack.key());
592 
593  // use recob::Track functionality as much as possible
594  // const double setlength = 2.5; double length = 0.0; // cm
595 
596  TVector3 pos_p = recotrack->LocationAtPoint<TVector3>(0);
597  float px = pos_p.X();
598  if (px > 0) px -= t0Corr(evt,detProp, pvtx);
599  else px += t0Corr(evt, detProp, pvtx);
600  pos_p.SetX(px);
601 
602  lep_dist = std::sqrt(pma::Dist2(primary, pos_p));
603  lep_distx = primary.X() - pos_p.X();
604  lep_disty = primary.Y() - pos_p.Y();
605  lep_distz = primary.Z() - pos_p.Z();
606  if (mclist.size())
607  {
608  lep_distreco = std::sqrt(pma::Dist2(closestvtx, pos_p));
609  lep_distrecox = closestvtx.X() - pos_p.X();
610  lep_distrecoy = closestvtx.Y() - pos_p.Y();
611  lep_distrecoz = closestvtx.Z() - pos_p.Z();
612  }
613 
614  cos = recoini.GetCos();
615  coslrx = mcdir3d.X(); coslry = mcdir3d.Y(); coslrz = mcdir3d.Z();
616 
617  /*************************************************************/
618  /* WARNING */
619  /*************************************************************/
620  /* The dQdx information in recob::Track has been deprecated */
621  /* since 2016 and in 11/2018 the recob::Track interface was */
622  /* changed and DQdxAtPoint and NumberdQdx were removed. */
623  /* Therefore the code below is now commented out */
624  /* (note that it was most likely not functional anyways). */
625  /* For any issue please contact: larsoft-team@fnal.gov */
626  /*************************************************************/
627  /*
628  size_t fp = 0; bool hitcoll = false;
629  for (size_t p = 0; p < recotrack->NumberTrajectoryPoints(); ++p)
630  if (recotrack->DQdxAtPoint(p, geo::kZ) > 0) {pos_p = recotrack->LocationAtPoint<TVector3>(p); fp = p; hitcoll = true; break;}
631 
632  // loop over trajectory point to get dQdx.
633  if (hitcoll)
634  for (size_t p = (fp+1); p < recotrack->NumberTrajectoryPoints(); ++p)
635  {
636  TVector3 pos = recotrack->LocationAtPoint<TVector3>(p);
637  length += std::sqrt(pma::Dist2(pos_p, pos));
638  pos_p = recotrack->LocationAtPoint<TVector3>(p);
639 
640  if (length > setlength) break;
641  dx = length;
642  double dqdx_p = recotrack->DQdxAtPoint(p, geo::kZ);
643  if (dqdx_p > 0) lep_dedx += recoinihit[p]->SummedADC();
644  }
645  if (dx > 0.0) lep_dedx /= dx;
646  */
647  /*************************************************************/
648  }
649 
650  }
651  }
652  }
653 }
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:219
int PdgCode() const
Definition: MCParticle.h:212
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:37
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
Point_t const & LocationAtPoint(size_t i) const
Definition: Track.h:126
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
T abs(T value)
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
std::string fGenieGenModuleLabel
bool insideFidVol(TLorentzVector const &pvtx) const
key_type key() const noexcept
Definition: Ptr.h:216
TVector3 findElDir(art::Ptr< simb::MCTruth > const mctruth) const
float t0Corr(art::Event const &evt, detinfo::DetectorPropertiesData const &detProp, TLorentzVector const &pvtx)
std::string fHitsModuleLabel
std::string fVertexModuleLabel
std::string fTrackModuleLabel
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
std::vector< TVector3 > findPhDir() const

Member Data Documentation

Float_t dunefd::IniSegReco::convdist
private

Definition at line 135 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::cos
private

Definition at line 122 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::coslrx
private

Definition at line 124 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::coslry
private

Definition at line 124 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::coslrz
private

Definition at line 124 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::coslx
private

Definition at line 123 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::cosly
private

Definition at line 123 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::coslz
private

Definition at line 123 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::dx
private

Definition at line 113 of file IniSegReco_module.cc.

int dunefd::IniSegReco::event
private

Definition at line 108 of file IniSegReco_module.cc.

std::string dunefd::IniSegReco::fClusterModuleLabel
private

Definition at line 138 of file IniSegReco_module.cc.

double dunefd::IniSegReco::fFidVolCut
private

Definition at line 143 of file IniSegReco_module.cc.

std::string dunefd::IniSegReco::fGenieGenModuleLabel
private

Definition at line 141 of file IniSegReco_module.cc.

std::string dunefd::IniSegReco::fHitsModuleLabel
private

Definition at line 137 of file IniSegReco_module.cc.

std::ofstream dunefd::IniSegReco::file
private

Definition at line 146 of file IniSegReco_module.cc.

std::ofstream dunefd::IniSegReco::file1
private

Definition at line 147 of file IniSegReco_module.cc.

pma::ProjectionMatchingAlg dunefd::IniSegReco::fProjectionMatchingAlg
private

Definition at line 144 of file IniSegReco_module.cc.

std::string dunefd::IniSegReco::fTrackModuleLabel
private

Definition at line 140 of file IniSegReco_module.cc.

TTree* dunefd::IniSegReco::fTree
private

Definition at line 104 of file IniSegReco_module.cc.

std::string dunefd::IniSegReco::fVertexModuleLabel
private

Definition at line 139 of file IniSegReco_module.cc.

short dunefd::IniSegReco::isdata
private

Definition at line 109 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::lep_dedx
private

Definition at line 112 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::lep_dist
private

Definition at line 114 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::lep_distreco
private

Definition at line 118 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::lep_distrecox
private

Definition at line 119 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::lep_distrecoy
private

Definition at line 120 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::lep_distrecoz
private

Definition at line 121 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::lep_distx
private

Definition at line 115 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::lep_disty
private

Definition at line 116 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::lep_distz
private

Definition at line 117 of file IniSegReco_module.cc.

int dunefd::IniSegReco::ngamma
private

Definition at line 134 of file IniSegReco_module.cc.

std::vector< pma::Track3D* > dunefd::IniSegReco::pmatracks
private

Definition at line 102 of file IniSegReco_module.cc.

int dunefd::IniSegReco::run
private

Definition at line 106 of file IniSegReco_module.cc.

int dunefd::IniSegReco::subrun
private

Definition at line 107 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::t0
private

Definition at line 125 of file IniSegReco_module.cc.

int dunefd::IniSegReco::trkindex
private

Definition at line 110 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::vtxrecomc
private

Definition at line 126 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::vtxrecomcx
private

Definition at line 127 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::vtxrecomcy
private

Definition at line 128 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::vtxrecomcz
private

Definition at line 129 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::vtxupstream
private

Definition at line 130 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::vtxupstreamx
private

Definition at line 131 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::vtxupstreamy
private

Definition at line 132 of file IniSegReco_module.cc.

Float_t dunefd::IniSegReco::vtxupstreamz
private

Definition at line 133 of file IniSegReco_module.cc.


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