13 #include "art_root_io/TFileService.h" 35 #include "TLorentzVector.h" 72 double fMinx{std::numeric_limits<double>::max()};
73 double fMaxx{std::numeric_limits<double>::min()};
74 double fMiny{std::numeric_limits<double>::max()};
75 double fMaxy{std::numeric_limits<double>::min()};
76 double fMinz{std::numeric_limits<double>::max()};
77 double fMaxz{std::numeric_limits<double>::min()};
112 fMaxx = std::max(fMaxx, tpcg.MaxX());
113 fMiny = std::min(fMiny, tpcg.MinY());
114 fMaxy = std::max(fMaxy, tpcg.MaxY());
115 fMinz = std::min(fMinz, tpcg.MinZ());
116 fMaxz = std::max(fMaxz, tpcg.MaxZ());
142 if (particle->
Process() !=
"primary")
continue;
144 TLorentzVector posvec = particle->
Position();
145 TVector3 pose(posvec.X(), posvec.Y(), posvec.Z());
148 if (particle->
PdgCode() == 111) {
151 TLorentzVector posvec3 = particle->
Position();
152 TVector3 pospi0(posvec3.X(), posvec3.Y(), posvec3.Z());
158 if (daughter1->
PdgCode() != 22)
continue;
161 if (daughter2->
PdgCode() != 22)
continue;
171 TVector3 mom1vec3(mom1.Px(), mom1.Py(), mom1.Pz());
173 TVector3 mom2vec3(mom2.Px(), mom2.Py(), mom2.Pz());
185 TVector3 vecnorm1 = mom1vec3.Unit();
187 TVector3 vecnorm2 = mom2vec3.Unit();
203 double dista = fabs(
fMinx - pvtx.X());
204 double distb = fabs(pvtx.X() -
fMaxx);
208 dista = fabs(
fMaxy - pvtx.Y());
209 distb = fabs(pvtx.Y() -
fMiny);
217 dista = fabs(
fMaxz - pvtx.Z());
218 distb = fabs(pvtx.Z() -
fMinz);
239 void endJob()
override;
246 TVector3
const& convmc,
332 fEvTree = tfs->make<TTree>(
"MultiShowers",
"showers3d");
350 fRecoTree = tfs->make<TTree>(
"Cascades",
"conv points");
362 fShTree = tfs->make<TTree>(
"Shower",
"conv point");
375 log <<
"******************** fEvFidVol = " <<
fEvFidVol <<
"\n";
376 log <<
"******************** fEvGMomCut = " <<
fEvGMomCut <<
"\n";
377 log <<
"******************** fEvComp = " <<
fEvComp <<
"\n";
378 log <<
"******************** fEvReco = " <<
fEvReco <<
"\n";
379 log <<
"******************** fEvInput = " <<
fEvInput <<
"\n";
380 log <<
"******************** fEv2Groups = " <<
fEv2Groups <<
"\n";
381 log <<
"******************** fEv2Good = " <<
fEv2Good <<
"\n";
383 log <<
"******************** reco % = " << double(fEvReco) / double(fEvInput) <<
"\n";
425 fMcth = 180.0F * (std::acos(cosinemc)) / TMath::Pi();
429 const double maxdist = 2.0;
462 for (
size_t s = 0; s < shsListHandle->size(); ++s) {
464 double mindist = maxdist;
467 for (
int i = 0; i <
fNgammas; ++i) {
469 if ((dist < mindist) && (idph != i)) {
487 for (
size_t s = 0; s < shsListHandle->size(); ++s) {
493 std::vector<double>
const& vecdedx = sh.
dEdx();
495 if (vecdedx.size() == 3) {
506 for (
size_t s = 0; s < shsListHandle->size(); ++s) {
508 double mindist = maxdist;
511 if (dist < mindist) {
516 std::vector<double> vecdedx = sh.
dEdx();
517 if (vecdedx.size() == 3) {
539 std::vector<std::pair<TVector3, TVector3>> lines;
543 std::pair<TVector3, TVector3> frontback1(sh1.
ShowerStart(),
545 std::pair<TVector3, TVector3> frontback2(sh2.
ShowerStart(),
547 lines.push_back(frontback1);
548 lines.push_back(frontback2);
557 if ((dist1_0 > dist2_0) || (dist1_1 > dist2_1)) {}
564 fRecth = 180.0F * (std::acos(cosine_reco)) / TMath::Pi();
597 geo::Point_t const vtx{convp[0].X(), convp[0].Y(), convp[0].Z()};
607 std::map<size_t, std::vector<size_t>> used;
613 double maxdist = 1.0;
614 if (geom->
HasTPC(idtpc)) {
619 for (
auto const& plane : readoutGeom.Iterate<
geo::PlaneGeo>(cryoid)) {
621 double mindist = maxdist;
623 for (
size_t c = 0; c < cluListHandle->size(); ++c) {
626 for (
auto const& ids : used)
627 for (
auto i : ids.second)
628 if (i == c) exist =
true;
631 std::vector<art::Ptr<recob::Hit>>
hits = fbc.at(c);
632 if (hits.size() < 20)
continue;
633 if (hits[0]->
WireID().Plane != plane.View())
continue;
637 if (dist < mindist) {
642 if (mindist < maxdist) used[conv].push_back(clid);
651 for (
auto const& ids : used) {
652 if (ids.second.size() > 1)
665 TVector3
const& convmc,
670 double mindist = 9999;
675 for (
size_t h = 0; h < v.size(); ++h) {
676 if ((v[h]->
WireID().
Plane == view) && (v[h]->WireID().TPC == tpc)) {
681 if (dist < mindist) { mindist =
dist; }
const TVector3 & ShowerStart() const
TVector3 const & GetDirgamma2() const
const TLorentzVector & Position(const int i=0) const
double SolveLeastSquares3D(const std::vector< std::pair< TVector3, TVector3 >> &lines, TVector3 &result)
double Dist2(const TVector2 &v1, const TVector2 &v2)
const simb::MCParticle * TrackIdToParticle_P(int id) const
MultiEMShowers(fhicl::ParameterSet const &p)
Declaration of signal hit object.
list_type::const_iterator const_iterator
art::InputTag fTrk3DModuleLabel
Geometry information for a single TPC.
void analyze(art::Event const &e) override
bool convCluster(art::Event const &evt)
CryostatID_t Cryostat
Index of cryostat.
TVector2 WireDriftToCm(detinfo::DetectorPropertiesData const &detProp, unsigned int wire, float drift, unsigned int plane, unsigned int tpc, unsigned int cryo)
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
std::string Process() const
double GetMomGamma2() const
int NumberDaughters() const
int Daughter(const int i) const
double getMinDist(detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &v, TVector3 const &convmc, size_t view, size_t tpc, size_t cryo)
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
art::InputTag fShsModuleLabel
art::InputTag fVtxModuleLabel
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
const std::vector< double > & dEdx() const
#define DEFINE_ART_MODULE(klass)
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
std::string EndProcess() const
art::InputTag fHitsModuleLabel
TVector3 const & GetPosgamma1() const
double P(const int i=0) const
T get(std::string const &key) const
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
const TVector3 & Direction() const
The data type to uniquely identify a TPC.
Declaration of cluster object.
TVector3 const & GetPrimary() const
TVector3 const & GetDirgamma1() const
const sim::ParticleList & ParticleList() const
CryostatID PositionToCryostatID(Point_t const &point) const
Returns the ID of the cryostat at specified location.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Implementation of the Projection Matching Algorithm.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
TVector3 const & GetPospi0() const
double GetMomGamma1() const
double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Provides recob::Track data product.
TVector2 GetProjectionToPlane(const TVector3 &p, unsigned int plane, unsigned int tpc, unsigned int cryo)
art::InputTag fCluModuleLabel
EventNumber_t event() const
TPCID_t TPC
Index of the TPC within its cryostat.
TVector3 const & GetPosgamma2() const
recob::tracking::Plane Plane
bool HasTPC(TPCID const &tpcid) const
Returns whether we have the specified TPC.
art framework interface to geometry description
The data type to uniquely identify a cryostat.
bool insideFidVol(const TLorentzVector &pvtx) const