238 mf::LogError(
"ClusterAna")<<
"Too many planes for this code...";
245 std::vector< art::Ptr<recob::Hit> > allhits;
247 if(allhits.size() == 0)
return;
253 if(clusterListHandle->size() == 0) {
254 std::cout<<
"No clusters in event. Hits size "<<allhits.size()<<
"\n";
261 double xyz[3] = {0,0,0};
263 if(vertexListHandle.
isValid()) {
265 for(
unsigned int ii = 0; ii < vertexListHandle->size(); ++ii){
277 std::vector<const simb::MCParticle*> plist2;
279 std::vector<std::vector<art::Ptr<recob::Hit>>> hlist2;
281 std::vector<std::vector<short>> truToCl;
283 std::vector<std::vector<unsigned short>> nTruHitInCl;
285 std::vector<unsigned short> nRecHitInCl;
290 float aveNuEP2mu = 0.;
291 float numNuEP2mu = 0.;
292 float aveNuEP2nm = 0.;
293 float numNuEP2nm = 0.;
299 int neutTrackID = -1;
300 std::vector<int> tidlist;
301 float neutEnergy = -1.;
302 int neutIntType = -1;
304 bool isNeutrino =
false;
305 bool isSingleParticle =
false;
309 int pdg = abs(part->
PdgCode());
315 float KE = 1000 * (part->
E() - part->
Mass());
319 if((isNeutrino || isSingleParticle) && neutTrackID < 0) {
320 neutTrackID = trackID;
323 neutEnergy = 1000. * theNeutrino.
Nu().
E();
325 neutCCNC = theNeutrino.
CCNC();
329 for(
unsigned short iv = 0; iv < recoVtxList.
size(); ++iv) {
330 recoVtxList[iv]->XYZ(xyz);
338 <<
"Origin "<<theTruth->
Origin()<<
" trackID "<<trackID
340 <<
" KE "<<(int)KE<<
" MeV " 341 <<
" neutTrackID "<<neutTrackID
342 <<
" Mother "<<part->
Mother()
345 bool isCharged = (pdg == 11) || (pdg == 13) || (pdg == 211)
346 || (pdg == 321) || (pdg == 2212);
349 if(!isCharged)
continue;
355 if(part->
Process() !=
"primary")
continue;
375 if(part->
Process() ==
"NeutronInelastic")
continue;
376 plist2.push_back(part);
377 tidlist.push_back(trackID);
379 std::vector<short> temp {-1, -1, -1};
380 truToCl.push_back(temp);
382 std::vector<unsigned short> temp2(3);
383 nTruHitInCl.push_back(temp2);
386 <<
"plist2["<<plist2.size() - 1<<
"]" 387 <<
" Origin "<<theTruth->
Origin()<<
" trackID "<<trackID
390 <<
" Mother "<<part->
Mother() + neutTrackID - 1
394 if(plist2.size() == 0)
return;
398 if(hlist2.size() != plist2.size()) {
400 <<
"MC particle list size "<<plist2.size()
401 <<
" != size of MC particle true hits lists "<<hlist2.size();
407 std::vector<std::pair<unsigned short, unsigned short>> moda;
412 unsigned short ii, dpl, jj, jpl, kpl;
413 for(ii = 0; ii < plist2.size(); ++ii) {
414 dpl = plist2.size() - 1 - ii;
416 if(plist2[dpl]->Mother() == 0)
continue;
418 if(abs(plist2[dpl]->PdgCode()) == 11)
continue;
420 int motherID = neutTrackID + plist2[dpl]->Mother() - 1;
422 if(motherID != neutTrackID)
continue;
425 for(kpl = 0; kpl < plist2.size(); ++kpl) {
426 if(plist2[kpl]->Mother() == motherID) ++ndtr;
429 if(ndtr > 1)
continue;
432 for(jj = 0; jj < plist2.size(); ++jj) {
433 jpl = plist2.size() - 1 - jj;
434 if(plist2[jpl]->TrackId() == motherID) {
441 mf::LogVerbatim(
"ClusterAna")<<
" mother of daughter "<<dpl<<
" not found. mpl = "<<mpl;
445 if(plist2[dpl]->PdgCode() != plist2[mpl]->PdgCode())
continue;
446 moda.push_back(std::make_pair(mpl, dpl));
452 for (
unsigned int ii = 0; ii < clusterListHandle->size(); ++ii){
458 nRecHitInCl.resize(clusters.
size());
462 std::map< geo::View_t, unsigned int > ViewToPlane;
463 for(
unsigned int plane=0; plane < geom->
Nplanes(); ++plane){
465 ViewToPlane[view] = plane;
467 for(
size_t icl = 0; icl < clusters.
size(); ++icl){
468 unsigned int plane = ViewToPlane[clusters[icl]->View()];
469 std::vector< art::Ptr<recob::Hit> > cluhits = fmh.at(icl);
471 nRecHitInCl[icl] = cluhits.size();
473 std::vector<unsigned short> nHitInPl2(plist2.size());
474 for(
size_t iht = 0; iht < cluhits.size(); ++iht){
480 for(
unsigned short ipl = 0; ipl < plist2.size(); ++ipl) {
481 unsigned short imat = 0;
482 for(imat = 0; imat < hlist2[ipl].size(); ++imat) {
483 if(cluhits[iht] == hlist2[ipl][imat])
break;
485 if(imat < hlist2[ipl].size()) {
490 if(hitInPl2 < 0)
continue;
494 for(
unsigned short imd = 0; imd < moda.size(); ++imd) {
495 if(moda[imd].second == hitInPl2) hitInPl2 = moda[imd].first;
498 ++nHitInPl2[hitInPl2];
502 unsigned short nhit = 0;
504 for(
unsigned int ipl = 0; ipl < nHitInPl2.size(); ++ipl) {
505 if(nHitInPl2[ipl] > nhit) {
506 nhit = nHitInPl2[ipl];
514 if(nhit > nTruHitInCl[imtru][plane]) {
515 truToCl[imtru][plane] = icl;
516 nTruHitInCl[imtru][plane] = nhit;
522 for(
unsigned short ipl = 0; ipl < plist2.size(); ++ipl) {
525 for(
unsigned short ii = 0; ii < moda.size(); ++ii) {
526 if(moda[ii].second == ipl) {
534 if(hlist2[ipl].size() < 3)
continue;
536 int trackID = plist2[ipl]->TrackId();
539 float KE = 1000 * (plist2[ipl]->E() - plist2[ipl]->Mass());
540 int PDG = abs(plist2[ipl]->PdgCode());
542 std::vector<short> nTru(geom->
Nplanes());
543 std::vector<short> nRec(geom->
Nplanes());
544 std::vector<short> nTruRec(geom->
Nplanes());
545 std::vector<float> eff(geom->
Nplanes());
546 std::vector<float> pur(geom->
Nplanes());
547 std::vector<float> ep(geom->
Nplanes());
548 for(
unsigned int plane = 0; plane < geom->
Nplanes(); ++plane) {
551 for(
unsigned short ii = 0; ii < hlist2[ipl].size(); ++ii){
552 if(ViewToPlane[hlist2[ipl][ii]->View()] == plane) ++nTru[plane];
556 unsigned short mom = ipl;
557 std::vector<std::pair<unsigned short, unsigned short>>::reverse_iterator
559 while(rit != moda.rend()) {
560 if((*rit).first == mom) {
561 unsigned short dau = (*rit).second;
562 for(
unsigned short jj = 0; jj < hlist2[dau].size(); ++jj) {
563 if(ViewToPlane[hlist2[dau][jj]->View()] == plane) ++nTru[plane];
574 if(nTru[plane] == 0)
continue;
576 if(nTru[plane] < 3)
continue;
577 short icl = truToCl[ipl][plane];
581 nRec[plane] = nRecHitInCl[icl];
583 nTruRec[plane] = nTruHitInCl[ipl][plane];
586 if(nTru[plane] > 0) eff[plane] = (float)nTruRec[plane] / (
float)nTru[plane];
587 if(nRec[plane] > 0) pur[plane] = (float)nTruRec[plane] / (
float)nRec[plane];
589 if(eff[plane] == 0) { eff[plane] = 0.001; }
590 if(pur[plane] == 0) { pur[plane] = 0.001; }
591 if(eff[plane] >= 1) { eff[plane] = 0.999; }
592 if(pur[plane] >= 1) { pur[plane] = 0.999; }
593 ep[plane] = eff[plane] * pur[plane];
594 if(
fPrintLevel == -1)
outFile<<
moduleID<<
" "<<evt.
id().
event()<<
" "<<trackID<<
" "<<PDG<<
" "<<(int)KE<<
" "<<plane<<
" "<<nTru[plane]<<
" "<<std::setprecision(3)<<eff[plane]<<
" "<<pur[plane]<<
"\n";
598 std::vector<float> temp;
600 std::sort(temp.begin(), temp.end());
602 unsigned short ii = temp.size() - 2;
603 float ep2 = temp[ii];
606 short ep2Cluster = 0;
607 for(
unsigned short jj = 0; jj < temp.size(); ++jj) {
610 ep2Cluster = truToCl[ipl][ep2Plane];
615 std::array<double, 2> clBeg, clEnd;
616 if(ep2Cluster >= 0) {
617 clBeg[0] = clusters[ep2Cluster]->StartWire();
618 clBeg[1] = clusters[ep2Cluster]->StartTick();
619 clEnd[0] = clusters[ep2Cluster]->EndWire();
620 clEnd[1] = clusters[ep2Cluster]->EndTick();
629 fCRE->Fill(eff[ep2Plane]);
630 fCRP->Fill(pur[ep2Plane]);
634 <<
">>>CREP2 "<<std::fixed<<std::setprecision(2)<<ep2
635 <<
" E "<<eff[ep2Plane]<<std::setprecision(2)<<
" P "<<pur[ep2Plane]
636 <<
" P:W:T "<<ep2Plane<<
":"<<(int)clBeg[0]<<
":"<<(
int)clBeg[1]
637 <<"-"<<ep2Plane<<":"<<(
int)clEnd[0]<<":"<<(
int)clEnd[1]
638 <<" PDG "<<PDG<<" KE "<<(
int)KE<<" MeV";
645 aveNuEP2mu += ep2 * wght;
648 aveNuEP2nm += ep2 * wght;
657 }
else if(PDG == 13) {
663 }
else if(PDG == 211) {
669 }
else if(PDG == 321) {
675 }
else if(PDG == 2212) {
683 <<
">>>NuEP2 "<<std::fixed<<std::setprecision(2)<<ep2
684 <<
" E "<<eff[ep2Plane]<<std::setprecision(2)<<
" P "<<pur[ep2Plane]
685 <<
" P:W:T "<<ep2Plane<<
":"<<(int)clBeg[0]<<
":"<<(
int)clBeg[1]
686 <<"-"<<ep2Plane<<":"<<(
int)clEnd[0]<<":"<<(
int)clEnd[1]
687 <<" PDG "<<PDG<<" KE "<<(
int)KE<<" MeV ";
691 mfp<<
" Truth P:W:T ";
692 for(
unsigned int plane = 0; plane < geom->
Nplanes(); ++plane) {
693 unsigned short loW = 9999;
694 unsigned short loT = 0;
695 unsigned short hiW = 0;
696 unsigned short hiT = 0;
697 for(
unsigned short ii = 0; ii < hlist2[ipl].size(); ++ii){
698 if(ViewToPlane[hlist2[ipl][ii]->View()] == plane) {
710 mfp<<plane<<
":"<<loW<<
":"<<loT<<
"-"<<plane<<
":"<<hiW<<
":"<<hiT<<
" ";
717 if(numNuEP2mu > 0.) ave1 = aveNuEP2mu/numNuEP2mu;
720 if(numNuEP2nm > 0.) ave2 = aveNuEP2nm/numNuEP2nm;
723 if(numCREP2 > 0.) ave3 = aveCREP2/numCREP2;
727 if(fPrintLevel > 0) {
728 std::string nuType =
"Other";
730 if(neutIntType == 1001) nuType =
"CCQE";
731 if(neutIntType == 1091) nuType =
"DIS";
732 if(neutIntType == 1097) nuType =
"COH";
733 if(neutIntType > 1002 && neutIntType < 1091) nuType =
"RES";
742 <<
" Enu "<<std::fixed<<std::setprecision(0)<<neutEnergy
743 <<
std::right<<std::fixed<<std::setprecision(2)
746 <<
" CosmicRays "<<ave3
747 <<
" CCNC "<<neutCCNC<<
" IntType "<<neutIntType;
double E(const int i=0) const
const std::vector< std::vector< art::Ptr< recob::Hit > > > TrackIdsToHits_Ps(std::vector< int > const &tkIds, std::vector< art::Ptr< recob::Hit > > const &hitsIn)
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::string fHitsModuleLabel
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)
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
geo::WireID WireID() const
Initial tdc tick for hit.
const simb::MCParticle & Nu() const
list_type::const_iterator const_iterator
simb::Origin_t Origin() const
WireID_t Wire
Index of the wire within its plane.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
std::string Process() const
std::vector< float > fElecKERange
TProfile * fNuEP2_KE_muon
int InteractionType() const
View_t View() const
Which coordinate does this plane measure.
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
void push_back(Ptr< U > const &p)
single particles thrown at the detector
TProfile * fNuEP2_KE_kaon
std::vector< float > fProtKERange
std::vector< float > fKaonKERange
std::string fVertexModuleLabel
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 &productInstanceName, Handle< PROD > &result) const
TProfile * fNuEP2_KE_pion
double Vz(const int i=0) const
EventNumber_t event() const
const sim::ParticleList & ParticleList()
std::string fClusterModuleLabel
TProfile * fNuEP2_KE_elec
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Event generator information.
std::vector< float > fPionKERange
std::vector< float > fMuonKERange
double Vy(const int i=0) const
const art::Ptr< simb::MCTruth > & TrackIdToMCTruth_P(int const &id)
TProfile * fNuEP2_KE_prot