180 if (geom->
Nplanes() > 3)
return;
185 std::vector<art::Ptr<recob::Hit>> allhits;
187 if (
empty(allhits))
return;
193 if (clusterListHandle->size() == 0)
return;
199 double xyz[3] = {0, 0, 0};
200 for (
unsigned int ii = 0; ii < vertexListHandle->size(); ++ii) {
210 for (
unsigned int ii = 0; ii < PFPListHandle->size(); ++ii) {
221 std::vector<const simb::MCParticle*> plist2;
223 std::vector<std::vector<art::Ptr<recob::Hit>>> hlist2;
225 std::vector<std::vector<short>> truToCl;
227 std::vector<std::vector<unsigned short>> nTruHitInCl;
229 std::vector<unsigned short> nRecHitInCl;
234 float aveNuEP2mu = 0.;
235 float numNuEP2mu = 0.;
236 float aveNuEP2nm = 0.;
237 float numNuEP2nm = 0.;
243 int neutTrackID = -1;
244 std::vector<int> tidlist;
245 float neutEnergy = -1.;
246 int neutIntType = -1;
261 << trackID <<
" PDG " << part->
PdgCode() <<
" E " << part->
E()
262 <<
" mass " << part->
Mass() <<
" Mother " 263 << part->
Mother() + neutTrackID <<
" Proc " << part->
Process();
268 neutTrackID = trackID;
270 neutEnergy = 1000. * theNeutrino.
Nu().
E();
272 neutCCNC = theNeutrino.
CCNC();
273 for (
unsigned short iv = 0; iv < recoVtxList.
size(); ++iv) {
274 recoVtxList[iv]->XYZ(xyz);
281 bool isCharged = (pdg == 11) || (pdg == 13) || (pdg == 211) || (pdg == 321) || (pdg == 2212);
283 if (!isCharged)
continue;
285 float KE = 1000 * (part->
E() - part->
Mass());
290 if (part->
Process() !=
"primary")
continue;
310 if (part->
Process() ==
"NeutronInelastic")
continue;
311 plist2.push_back(part);
312 tidlist.push_back(trackID);
314 std::vector<short> temp{-1, -1, -1};
315 truToCl.push_back(temp);
317 std::vector<unsigned short> temp2(3);
318 nTruHitInCl.push_back(temp2);
322 <<
" trackID " << trackID <<
" PDG " << part->
PdgCode() <<
" KE " 323 << (int)KE <<
" Mother " << part->
Mother() + neutTrackID
324 <<
" Proc " << part->
Process();
327 if (
empty(plist2))
return;
331 if (hlist2.size() != plist2.size()) {
332 mf::LogError(
"PFPAna") <<
"MC particle list size " << plist2.size()
333 <<
" != size of MC particle true hits lists " << hlist2.size();
339 std::vector<std::pair<unsigned short, unsigned short>> moda;
344 for (
unsigned short dpl = plist2.size() - 1; dpl > 0; --dpl) {
346 if (plist2[dpl]->Mother() == 0)
continue;
348 if (
abs(plist2[dpl]->PdgCode()) == 11)
continue;
350 int motherID = neutTrackID + plist2[dpl]->Mother() - 1;
352 if (motherID < 0)
continue;
355 for (
unsigned short kpl = 0; kpl < plist2.size(); ++kpl) {
356 if (plist2[kpl]->Mother() == motherID) ++ndtr;
359 if (ndtr > 1)
continue;
362 for (
unsigned short jpl = dpl - 1; jpl > 0; --jpl) {
363 if (plist2[jpl]->TrackId() == motherID) {
369 if (mpl < 0)
continue;
371 if (plist2[dpl]->PdgCode() != plist2[mpl]->PdgCode())
continue;
372 moda.push_back(std::make_pair(mpl, dpl));
378 for (
unsigned int ii = 0; ii < clusterListHandle->size(); ++ii) {
384 nRecHitInCl.resize(clusters.
size());
388 std::map<geo::View_t, unsigned int> ViewToPlane;
391 ViewToPlane[view] = plane.ID().Plane;
393 for (
size_t icl = 0; icl < clusters.
size(); ++icl) {
394 unsigned int plane = ViewToPlane[clusters[icl]->View()];
395 std::vector<art::Ptr<recob::Hit>> cluhits = fmh.at(icl);
397 nRecHitInCl[icl] = cluhits.size();
399 std::vector<unsigned short> nHitInPl2(plist2.size());
400 for (
size_t iht = 0; iht < cluhits.size(); ++iht) {
403 for (
unsigned short ipl = 0; ipl < plist2.size(); ++ipl) {
404 unsigned short imat = 0;
405 for (imat = 0; imat < hlist2[ipl].size(); ++imat) {
406 if (cluhits[iht] == hlist2[ipl][imat])
break;
408 if (imat < hlist2[ipl].
size()) {
413 if (hitInPl2 < 0)
continue;
417 for (
unsigned short imd = 0; imd < moda.size(); ++imd) {
418 if (moda[imd].
second == hitInPl2) hitInPl2 = moda[imd].first;
421 ++nHitInPl2[hitInPl2];
425 unsigned short nhit = 0;
427 for (
unsigned int ipl = 0; ipl < nHitInPl2.size(); ++ipl) {
428 if (nHitInPl2[ipl] > nhit) {
429 nhit = nHitInPl2[ipl];
437 if (nhit > nTruHitInCl[imtru][plane]) {
438 truToCl[imtru][plane] = icl;
439 nTruHitInCl[imtru][plane] = nhit;
445 for (
unsigned short ipl = 0; ipl < plist2.size(); ++ipl) {
448 for (
unsigned short ii = 0; ii < moda.size(); ++ii) {
449 if (moda[ii].
second == ipl) {
454 if (skipit)
continue;
457 if (hlist2[ipl].
size() < 3)
continue;
459 int trackID = plist2[ipl]->TrackId();
462 float KE = 1000 * (plist2[ipl]->E() - plist2[ipl]->Mass());
463 int PDG =
abs(plist2[ipl]->PdgCode());
465 std::vector<short> nTru(geom->
Nplanes());
466 std::vector<short> nRec(geom->
Nplanes());
467 std::vector<short> nTruRec(geom->
Nplanes());
468 std::vector<float> eff(geom->
Nplanes());
469 std::vector<float> pur(geom->
Nplanes());
470 std::vector<float> ep(geom->
Nplanes());
471 for (
unsigned int plane = 0; plane < geom->
Nplanes(); ++plane) {
474 for (
unsigned short ii = 0; ii < hlist2[ipl].size(); ++ii) {
475 if (ViewToPlane[hlist2[ipl][ii]->View()] == plane) ++nTru[plane];
478 unsigned short mom = ipl;
479 std::vector<std::pair<unsigned short, unsigned short>>::reverse_iterator rit =
481 while (rit != moda.rend()) {
482 if ((*rit).first == mom) {
483 unsigned short dau = (*rit).second;
484 for (
unsigned short jj = 0; jj < hlist2[dau].size(); ++jj) {
485 if (ViewToPlane[hlist2[dau][jj]->View()] == plane) ++nTru[plane];
494 if (nTru[plane] == 0) {
continue; }
495 short icl = truToCl[ipl][plane];
496 nRec[plane] = nRecHitInCl[icl];
497 nTruRec[plane] = nTruHitInCl[ipl][plane];
498 if (nTru[plane] > 0) eff[plane] = (float)nTruRec[plane] / (
float)nTru[plane];
499 if (nRec[plane] > 0) pur[plane] = (float)nTruRec[plane] / (
float)nRec[plane];
500 ep[plane] = eff[plane] * pur[plane];
503 std::vector<float> temp;
505 std::sort(temp.begin(), temp.end());
507 unsigned short ii = temp.size() - 2;
508 float ep2 = temp[ii];
511 short ep2Cluster = 0;
512 for (
unsigned short jj = 0; jj < temp.size(); ++jj) {
515 ep2Cluster = truToCl[ipl][ep2Plane];
520 std::array<double, 2> clBeg, clEnd;
521 if (ep2Cluster >= 0) {
522 clBeg[0] = clusters[ep2Cluster]->StartWire();
523 clBeg[1] = clusters[ep2Cluster]->StartTick();
524 clEnd[0] = clusters[ep2Cluster]->EndWire();
525 clEnd[1] = clusters[ep2Cluster]->EndTick();
534 fCRE->Fill(eff[ep2Plane]);
535 fCRP->Fill(pur[ep2Plane]);
540 <<
">>>CREP2 " << std::fixed << std::setprecision(2) << ep2 <<
" E " << eff[ep2Plane]
541 << std::setprecision(2) <<
" P " << pur[ep2Plane] <<
" P:W:T " << ep2Plane <<
":" 542 << (int)clBeg[0] <<
":" << (
int)clBeg[1] << "-" << ep2Plane << ":" << (
int)clEnd[0]
543 << ":" << (
int)clEnd[1] << " PDG " << PDG << " KE " << (
int)KE << " MeV";
550 aveNuEP2mu += ep2 * wght;
554 aveNuEP2nm += ep2 * wght;
564 else if (PDG == 13) {
571 else if (PDG == 211) {
578 else if (PDG == 321) {
585 else if (PDG == 2212) {
594 <<
">>>NuEP2 " << std::fixed << std::setprecision(2) << ep2 <<
" E " << eff[ep2Plane]
595 << std::setprecision(2) <<
" P " << pur[ep2Plane] <<
" P:W:T " << ep2Plane <<
":" 596 << (int)clBeg[0] <<
":" << (
int)clBeg[1] << "-" << ep2Plane << ":" << (
int)clEnd[0]
597 << ":" << (
int)clEnd[1] << " PDG " << PDG << " KE " << (
int)KE << " MeV ";
601 mfp <<
" Truth P:W:T ";
602 for (
unsigned int plane = 0; plane < geom->
Nplanes(); ++plane) {
603 unsigned short loW = 9999;
604 unsigned short loT = 0;
605 unsigned short hiW = 0;
606 unsigned short hiT = 0;
607 for (
unsigned short ii = 0; ii < hlist2[ipl].size(); ++ii) {
608 if (ViewToPlane[hlist2[ipl][ii]->View()] == plane) {
620 mfp << plane <<
":" << loW <<
":" << loT <<
"-" << plane <<
":" << hiW <<
":" << hiT
628 if (numNuEP2mu > 0.) ave1 = aveNuEP2mu / numNuEP2mu;
631 if (numNuEP2nm > 0.) ave2 = aveNuEP2nm / numNuEP2nm;
634 if (numCREP2 > 0.) ave3 = aveCREP2 / numCREP2;
636 if (fPrintLevel > 0) {
637 std::string nuType =
"Other";
639 if (neutIntType == 1001) nuType =
"CCQE";
640 if (neutIntType == 1091) nuType =
"DIS";
641 if (neutIntType == 1097) nuType =
"COH";
642 if (neutIntType > 1002 && neutIntType < 1091) nuType =
"RES";
651 << std::fixed << std::setprecision(0) << neutEnergy <<
std::right 652 << std::fixed << std::setprecision(2) <<
" NuMuons " << ave1
653 <<
" NuPiKp " << ave2 <<
" CosmicRays " << ave3 <<
" CCNC " 654 << neutCCNC <<
" IntType " << neutIntType;
double E(const int i=0) const
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::vector< float > fElecKERange
const simb::MCNeutrino & GetNeutrino() const
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
std::string fVertexModuleLabel
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
const simb::MCParticle & Nu() const
list_type::const_iterator const_iterator
simb::Origin_t Origin() const
std::string fHitsModuleLabel
std::vector< float > fProtKERange
constexpr auto abs(T v)
Returns the absolute value of the argument.
WireID_t Wire
Index of the wire within its plane.
geo::WireID const & WireID() const
Initial tdc tick for hit.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
std::string Process() const
std::string fClusterModuleLabel
std::vector< float > fPionKERange
TProfile * fNuEP2_KE_kaon
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
int InteractionType() const
TProfile * fNuEP2_KE_elec
void push_back(Ptr< U > const &p)
const art::Ptr< simb::MCTruth > & TrackIdToMCTruth_P(int id) const
TProfile * fNuEP2_KE_pion
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
std::string fPFParticleModuleLabel
std::vector< std::vector< art::Ptr< recob::Hit > > > TrackIdsToHits_Ps(detinfo::DetectorClocksData const &clockData, std::vector< int > const &tkIds, std::vector< art::Ptr< recob::Hit >> const &hitsIn) const
The data type to uniquely identify a TPC.
const sim::ParticleList & ParticleList() const
TProfile * fNuEP2_KE_prot
double Vx(const int i=0) const
float PeakTime() const
Time of the signal peak, in tick units.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
TProfile * fNuEP2_KE_muon
double Vz(const int i=0) const
EventNumber_t event() const
std::vector< float > fMuonKERange
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
second_as<> second
Type of time stored in seconds, in double precision.
Event generator information.
double Vy(const int i=0) const
std::vector< float > fKaonKERange
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.