20 #include "art_root_io/TFileService.h" 136 if (found != std::string::npos)
moduleID = 1;
138 if (found != std::string::npos)
moduleID = 2;
140 if (found != std::string::npos)
moduleID = 3;
142 if (found != std::string::npos)
moduleID = 4;
157 fNClusters = tfs->make<TH1F>(
"fNoClustersInEvent",
"Number of Clusters", 40, 0, 400);
158 fNHitInCluster = tfs->make<TH1F>(
"fNHitInCluster",
"NHitInCluster", 100, 0, 100);
162 fCREP2 = tfs->make<TH1F>(
"CREP2",
"CREP2", 50, 0, 1);
163 fCRE = tfs->make<TH1F>(
"CRE",
"CR Efficiency", 50, 0, 1);
164 fCRP = tfs->make<TH1F>(
"CRP",
"CR Purity", 50, 0, 1);
167 fNuKE_elec = tfs->make<TH1F>(
"NuKE_elec",
"NuKE electron", 100, 0, 4000);
168 fNuKE_muon = tfs->make<TH1F>(
"NuKE_muon",
"NuKE muon", 100, 0, 4000);
169 fNuKE_pion = tfs->make<TH1F>(
"NuKE_pion",
"NuKE pion", 100, 0, 4000);
170 fNuKE_kaon = tfs->make<TH1F>(
"NuKE_kaon",
"NuKE kaon", 100, 0, 4000);
171 fNuKE_prot = tfs->make<TH1F>(
"NuKE_prot",
"NuKE proton", 100, 0, 4000);
173 fNuEP2_elec = tfs->make<TH1F>(
"NuEP2_elec",
"NuEP2 electron", 50, 0, 1);
174 fNuEP2_muon = tfs->make<TH1F>(
"NuEP2_muon",
"NuEP2 muon", 50, 0, 1);
175 fNuEP2_pion = tfs->make<TH1F>(
"NuEP2_pion",
"NuEP2 pion", 50, 0, 1);
176 fNuEP2_kaon = tfs->make<TH1F>(
"NuEP2_kaon",
"NuEP2 kaon", 50, 0, 1);
177 fNuEP2_prot = tfs->make<TH1F>(
"NuEP2_prot",
"NuEP2 proton", 50, 0, 1);
179 fNuE_elec = tfs->make<TH1F>(
"NuE_elec",
"Nu Efficiency electron", 50, 0, 1);
180 fNuE_muon = tfs->make<TH1F>(
"NuE_muon",
"Nu Efficiency muon", 50, 0, 1);
181 fNuE_pion = tfs->make<TH1F>(
"NuE_pion",
"Nu Efficiency pion", 50, 0, 1);
182 fNuE_kaon = tfs->make<TH1F>(
"NuE_kaon",
"Nu Efficiency kaon", 50, 0, 1);
183 fNuE_prot = tfs->make<TH1F>(
"NuE_prot",
"Nu Efficiency proton", 50, 0, 1);
185 fNuP_elec = tfs->make<TH1F>(
"NuP_elec",
"Nu Purity electron", 50, 0, 1);
186 fNuP_muon = tfs->make<TH1F>(
"NuP_muon",
"Nu Purity muon", 50, 0, 1);
187 fNuP_pion = tfs->make<TH1F>(
"NuP_pion",
"Nu Purity pion", 50, 0, 1);
188 fNuP_kaon = tfs->make<TH1F>(
"NuP_kaon",
"Nu Purity kaon", 50, 0, 1);
189 fNuP_prot = tfs->make<TH1F>(
"NuP_prot",
"Nu Purity proton", 50, 0, 1);
192 fNuVtx_dx = tfs->make<TH1F>(
"Vtx dx",
"Vtx dx", 80, -10, 10);
193 fNuVtx_dy = tfs->make<TH1F>(
"Vtx dy",
"Vtx dy", 80, -10, 10);
194 fNuVtx_dz = tfs->make<TH1F>(
"Vtx dz",
"Vtx dz", 80, -10, 10);
196 fNuEP2_KE_elec = tfs->make<TProfile>(
"NuEP2_KE_elec",
"NuEP2 electron vs KE", 200, 0, 2000);
197 fNuEP2_KE_muon = tfs->make<TProfile>(
"NuEP2_KE_muon",
"NuEP2 muon vs KE", 200, 0, 2000);
198 fNuEP2_KE_pion = tfs->make<TProfile>(
"NuEP2_KE_pion",
"NuEP2 pion vs KE", 200, 0, 2000);
199 fNuEP2_KE_kaon = tfs->make<TProfile>(
"NuEP2_KE_kaon",
"NuEP2 kaon vs KE", 200, 0, 2000);
200 fNuEP2_KE_prot = tfs->make<TProfile>(
"NuEP2_KE_prot",
"NuEP2 proton vs KE", 200, 0, 2000);
213 mf::LogError(
"ClusterAna") <<
"Too many planes for this code...";
220 std::vector<art::Ptr<recob::Hit>> allhits;
222 if (allhits.size() == 0)
return;
228 if (clusterListHandle->size() == 0) {
229 std::cout <<
"No clusters in event. Hits size " << allhits.size() <<
"\n";
235 double xyz[3] = {0, 0, 0};
237 if (vertexListHandle.
isValid()) {
239 for (
unsigned int ii = 0; ii < vertexListHandle->size(); ++ii) {
251 std::vector<const simb::MCParticle*> plist2;
253 std::vector<std::vector<art::Ptr<recob::Hit>>> hlist2;
255 std::vector<std::vector<short>> truToCl;
257 std::vector<std::vector<unsigned short>> nTruHitInCl;
259 std::vector<unsigned short> nRecHitInCl;
264 float aveNuEP2mu = 0.;
265 float numNuEP2mu = 0.;
266 float aveNuEP2nm = 0.;
267 float numNuEP2nm = 0.;
273 int neutTrackID = -1;
274 std::vector<int> tidlist;
275 float neutEnergy = -1.;
276 int neutIntType = -1;
278 bool isNeutrino =
false;
279 bool isSingleParticle =
false;
292 float KE = 1000 * (part->
E() - part->
Mass());
296 if ((isNeutrino || isSingleParticle) && neutTrackID < 0) {
297 neutTrackID = trackID;
300 neutEnergy = 1000. * theNeutrino.
Nu().
E();
302 neutCCNC = theNeutrino.
CCNC();
305 for (
unsigned short iv = 0; iv < recoVtxList.
size(); ++iv) {
306 recoVtxList[iv]->XYZ(xyz);
315 <<
" PDG " << part->
PdgCode() <<
" KE " << (int)KE <<
" MeV " 316 <<
" neutTrackID " << neutTrackID <<
" Mother " 319 bool isCharged = (pdg == 11) || (pdg == 13) || (pdg == 211) || (pdg == 321) || (pdg == 2212);
321 if (!isCharged)
continue;
327 if (part->
Process() !=
"primary")
continue;
347 if (part->
Process() ==
"NeutronInelastic")
continue;
348 plist2.push_back(part);
349 tidlist.push_back(trackID);
351 std::vector<short> temp{-1, -1, -1};
352 truToCl.push_back(temp);
354 std::vector<unsigned short> temp2(3);
355 nTruHitInCl.push_back(temp2);
359 <<
"plist2[" << plist2.size() - 1 <<
"]" 360 <<
" Origin " << theTruth->
Origin() <<
" trackID " << trackID <<
" PDG " 361 << part->
PdgCode() <<
" KE " << (int)KE <<
" Mother " << part->
Mother() + neutTrackID - 1
362 <<
" Proc " << part->
Process();
365 if (plist2.size() == 0)
return;
369 if (hlist2.size() != plist2.size()) {
370 mf::LogError(
"ClusterAna") <<
"MC particle list size " << plist2.size()
371 <<
" != size of MC particle true hits lists " << hlist2.size();
377 std::vector<std::pair<unsigned short, unsigned short>> moda;
382 unsigned short ii, dpl, jj, jpl, kpl;
383 for (ii = 0; ii < plist2.size(); ++ii) {
384 dpl = plist2.size() - 1 - ii;
386 if (plist2[dpl]->Mother() == 0)
continue;
388 if (
abs(plist2[dpl]->PdgCode()) == 11)
continue;
390 int motherID = neutTrackID + plist2[dpl]->Mother() - 1;
392 if (motherID != neutTrackID)
continue;
395 for (kpl = 0; kpl < plist2.size(); ++kpl) {
396 if (plist2[kpl]->Mother() == motherID) ++ndtr;
399 if (ndtr > 1)
continue;
402 for (jj = 0; jj < plist2.size(); ++jj) {
403 jpl = plist2.size() - 1 - jj;
404 if (plist2[jpl]->TrackId() == motherID) {
412 <<
" mother of daughter " << dpl <<
" not found. mpl = " << mpl;
416 if (plist2[dpl]->PdgCode() != plist2[mpl]->PdgCode())
continue;
417 moda.push_back(std::make_pair(mpl, dpl));
423 for (
unsigned int ii = 0; ii < clusterListHandle->size(); ++ii) {
429 nRecHitInCl.resize(clusters.
size());
433 std::map<geo::View_t, unsigned int> ViewToPlane;
436 ViewToPlane[view] = plane.ID().Plane;
438 for (
size_t icl = 0; icl < clusters.
size(); ++icl) {
439 unsigned int plane = ViewToPlane[clusters[icl]->View()];
440 std::vector<art::Ptr<recob::Hit>> cluhits = fmh.at(icl);
442 nRecHitInCl[icl] = cluhits.size();
444 std::vector<unsigned short> nHitInPl2(plist2.size());
445 for (
size_t iht = 0; iht < cluhits.size(); ++iht) {
449 for (
unsigned short ipl = 0; ipl < plist2.size(); ++ipl) {
450 unsigned short imat = 0;
451 for (imat = 0; imat < hlist2[ipl].size(); ++imat) {
452 if (cluhits[iht] == hlist2[ipl][imat])
break;
454 if (imat < hlist2[ipl].
size()) {
459 if (hitInPl2 < 0)
continue;
463 for (
unsigned short imd = 0; imd < moda.size(); ++imd) {
464 if (moda[imd].
second == hitInPl2) hitInPl2 = moda[imd].first;
467 ++nHitInPl2[hitInPl2];
471 unsigned short nhit = 0;
473 for (
unsigned int ipl = 0; ipl < nHitInPl2.size(); ++ipl) {
474 if (nHitInPl2[ipl] > nhit) {
475 nhit = nHitInPl2[ipl];
483 if (nhit > nTruHitInCl[imtru][plane]) {
484 truToCl[imtru][plane] = icl;
485 nTruHitInCl[imtru][plane] = nhit;
491 for (
unsigned short ipl = 0; ipl < plist2.size(); ++ipl) {
494 for (
unsigned short ii = 0; ii < moda.size(); ++ii) {
495 if (moda[ii].
second == ipl) {
500 if (skipit)
continue;
503 if (hlist2[ipl].
size() < 3)
continue;
505 int trackID = plist2[ipl]->TrackId();
508 float KE = 1000 * (plist2[ipl]->E() - plist2[ipl]->Mass());
509 int PDG =
abs(plist2[ipl]->PdgCode());
511 std::vector<short> nTru(geom->
Nplanes());
512 std::vector<short> nRec(geom->
Nplanes());
513 std::vector<short> nTruRec(geom->
Nplanes());
514 std::vector<float> eff(geom->
Nplanes());
515 std::vector<float> pur(geom->
Nplanes());
516 std::vector<float> ep(geom->
Nplanes());
517 for (
unsigned int plane = 0; plane < geom->
Nplanes(); ++plane) {
520 for (
unsigned short ii = 0; ii < hlist2[ipl].size(); ++ii) {
521 if (ViewToPlane[hlist2[ipl][ii]->View()] == plane) ++nTru[plane];
525 unsigned short mom = ipl;
526 std::vector<std::pair<unsigned short, unsigned short>>::reverse_iterator rit =
528 while (rit != moda.rend()) {
529 if ((*rit).first == mom) {
530 unsigned short dau = (*rit).second;
531 for (
unsigned short jj = 0; jj < hlist2[dau].size(); ++jj) {
532 if (ViewToPlane[hlist2[dau][jj]->View()] == plane) ++nTru[plane];
541 if (nTru[plane] == 0)
continue;
543 if (nTru[plane] < 3)
continue;
544 short icl = truToCl[ipl][plane];
545 if (icl < 0) { nRec[plane] = 0; }
547 nRec[plane] = nRecHitInCl[icl];
549 nTruRec[plane] = nTruHitInCl[ipl][plane];
552 if (nTru[plane] > 0) eff[plane] = (float)nTruRec[plane] / (
float)nTru[plane];
553 if (nRec[plane] > 0) pur[plane] = (float)nTruRec[plane] / (
float)nRec[plane];
555 if (eff[plane] == 0) { eff[plane] = 0.001; }
556 if (pur[plane] == 0) { pur[plane] = 0.001; }
557 if (eff[plane] >= 1) { eff[plane] = 0.999; }
558 if (pur[plane] >= 1) { pur[plane] = 0.999; }
559 ep[plane] = eff[plane] * pur[plane];
562 << (int)KE <<
" " << plane <<
" " << nTru[plane] <<
" " << std::setprecision(3)
563 << eff[plane] <<
" " << pur[plane] <<
"\n";
566 std::vector<float> temp;
568 std::sort(temp.begin(), temp.end());
570 unsigned short ii = temp.size() - 2;
571 float ep2 = temp[ii];
574 short ep2Cluster = 0;
575 for (
unsigned short jj = 0; jj < temp.size(); ++jj) {
578 ep2Cluster = truToCl[ipl][ep2Plane];
583 std::array<double, 2> clBeg, clEnd;
584 if (ep2Cluster >= 0) {
585 clBeg[0] = clusters[ep2Cluster]->StartWire();
586 clBeg[1] = clusters[ep2Cluster]->StartTick();
587 clEnd[0] = clusters[ep2Cluster]->EndWire();
588 clEnd[1] = clusters[ep2Cluster]->EndTick();
597 fCRE->Fill(eff[ep2Plane]);
598 fCRP->Fill(pur[ep2Plane]);
603 <<
">>>CREP2 " << std::fixed << std::setprecision(2) << ep2 <<
" E " << eff[ep2Plane]
604 << std::setprecision(2) <<
" P " << pur[ep2Plane] <<
" P:W:T " << ep2Plane <<
":" 605 << (int)clBeg[0] <<
":" << (
int)clBeg[1] <<
"-" << ep2Plane <<
":" << (int)clEnd[0]
606 <<
":" << (
int)clEnd[1] <<
" PDG " << PDG <<
" KE " << (int)KE <<
" MeV";
613 aveNuEP2mu += ep2 * wght;
617 aveNuEP2nm += ep2 * wght;
627 else if (PDG == 13) {
634 else if (PDG == 211) {
641 else if (PDG == 321) {
648 else if (PDG == 2212) {
657 <<
">>>NuEP2 " << std::fixed << std::setprecision(2) << ep2 <<
" E " << eff[ep2Plane]
658 << std::setprecision(2) <<
" P " << pur[ep2Plane] <<
" P:W:T " << ep2Plane <<
":" 659 << (int)clBeg[0] <<
":" << (
int)clBeg[1] <<
"-" << ep2Plane <<
":" << (int)clEnd[0]
660 <<
":" << (
int)clEnd[1] <<
" PDG " << PDG <<
" KE " << (int)KE <<
" MeV ";
664 mfp <<
" Truth P:W:T ";
665 for (
unsigned int plane = 0; plane < geom->
Nplanes(); ++plane) {
666 unsigned short loW = 9999;
667 unsigned short loT = 0;
668 unsigned short hiW = 0;
669 unsigned short hiT = 0;
670 for (
unsigned short ii = 0; ii < hlist2[ipl].size(); ++ii) {
671 if (ViewToPlane[hlist2[ipl][ii]->View()] == plane) {
683 mfp << plane <<
":" << loW <<
":" << loT <<
"-" << plane <<
":" << hiW <<
":" << hiT
691 if (numNuEP2mu > 0.) ave1 = aveNuEP2mu / numNuEP2mu;
694 if (numNuEP2nm > 0.) ave2 = aveNuEP2nm / numNuEP2nm;
697 if (numCREP2 > 0.) ave3 = aveCREP2 / numCREP2;
700 std::string nuType =
"Other";
702 if (neutIntType == 1001) nuType =
"CCQE";
703 if (neutIntType == 1091) nuType =
"DIS";
704 if (neutIntType == 1097) nuType =
"COH";
705 if (neutIntType > 1002 && neutIntType < 1091) nuType =
"RES";
714 <<
"EvtEP2 " << evt.
id().
event() <<
" NuType " << nuType <<
" Enu " << std::fixed
715 << std::setprecision(0) << neutEnergy <<
std::right << std::fixed << std::setprecision(2)
716 <<
" NuMuons " << ave1 <<
" NuPiKp " << ave2 <<
" CosmicRays " << ave3 <<
" CCNC " 717 << neutCCNC <<
" IntType " << neutIntType;
double E(const int i=0) const
void XYZ(double *xyz) const
Legacy method to access vertex position, preserved to avoid breaking code. Please try to use Vertex::...
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
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)
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Declaration of signal hit object.
const simb::MCParticle & Nu() const
list_type::const_iterator const_iterator
simb::Origin_t Origin() const
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
Cluster finding and building.
EDAnalyzer(fhicl::ParameterSet const &pset)
std::string Process() const
void analyze(const art::Event &evt) override
std::vector< float > fElecKERange
bool isValid() const noexcept
TProfile * fNuEP2_KE_muon
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
int InteractionType() const
ClusterAna(fhicl::ParameterSet const &pset)
#define DEFINE_ART_MODULE(klass)
void push_back(Ptr< U > const &p)
const art::Ptr< simb::MCTruth > & TrackIdToMCTruth_P(int id) const
single particles thrown at the detector
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
TProfile * fNuEP2_KE_kaon
std::vector< float > fProtKERange
std::vector< float > fKaonKERange
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
std::string fVertexModuleLabel
The data type to uniquely identify a TPC.
Declaration of cluster object.
const sim::ParticleList & ParticleList() const
double Vx(const int i=0) const
float PeakTime() const
Time of the signal peak, in tick units.
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Encapsulate the construction of a single detector plane.
TProfile * fNuEP2_KE_pion
double Vz(const int i=0) const
EventNumber_t event() const
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
std::string fClusterModuleLabel
TProfile * fNuEP2_KE_elec
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Particle list in DetSim contains Monte Carlo particle information.
second_as<> second
Type of time stored in seconds, in double precision.
Event generator information.
std::vector< float > fPionKERange
std::vector< float > fMuonKERange
double Vy(const int i=0) const
art framework interface to geometry description
TProfile * fNuEP2_KE_prot