37 #include "Utilities/AssociationUtil.h" 71 std::vector<cheat::NeutrinoEffPur>
sortEffPur(
const std::vector<cheat::NeutrinoEffPur>& nuEP);
304 fOutTree = tfs->
make<TTree>(
"SlicerAna",
"Slicer Analysis Tree");
416 "# of Nu;# of neutrinos in event;count",
419 "# of MCTruths;# of MCTruth objects in event;count",
422 "# of non-noise slices;# of non-noise slices per event;count",
425 "# of MCTruths that contributed to non-noise slices vs. # of non-noise slices;# of MCTruth objects;# of slices",
426 101,-0.5,100.5,101,-0.5,100.5);
428 "# of noise slices;# of noise slices per event;count",
431 "# of hits per non-noise slice;# of hits per non-noise slice;count",
434 "# of hits in the noise slice;# of hits per noise slice;count",
437 fComp = tfs->
make<TH1F>(
"SliceCompleteness",
438 "Completeness for most pure MCTruth in slice; slice completeness; count",
440 fPur = tfs->
make<TH1F>(
"SlicePurity",
441 "Purity for most pure MCTruth in slice; slice purity; count",
444 "Completeness/Purity for 'most pure' MCTruth by slice;Slice Purity;Slice Completeness",
445 101,-0.005,1.005,101,-0.005,1.005);
447 "Completeness for all MCTruths in NOISE slice; completeness in noise slice; count",
450 "Purity for all MCTruths in NOISE slice; purity in noise slice; count",
453 fNSvsVE =
new TH1F(
"SliceVisE",
"Slice VisE;slice visible energy (GeV)",80,0.0,20.0);
454 fNGvsVE =
new TH1F(
"GoodSliceVisE",
"good slice VisE;good slice (> 90% efficiency, > 90% purity) visible energy (GeV)",80,0.0,20.0);
458 "Efficiency vs. visible energy of best matched MCTruth for all non-noise slices;Visible energy [GeV];slice efficiency",
462 "Plane vs. X cell for all hits in non-noise slices;plane;X cell",
463 1320,-20.0,1300.0,510,-10.0,500.0);
465 "Plane vs. Y cell for all hits in non-noise slices;plane;Y cell",
466 1320,-20.0,1300.0,510,-10.0,500.0);
468 "Time for all hits in non-noise slices; time [ns]; count",
469 12000,-50000.0,550000.0);
471 "avg hit plane vs. avg hit X cell;avg slice plane;avg slice X cell",
472 660,-20.0,1300.0,255,-10.0,500.0);
474 "ave hit plane vs. ave hit Y cell;avg slice plane;avg slice Y cell",
475 660,-20.0,1300.0,255,-10.0,500.0);
477 "ave hit time per slice; avg slice time [ns]; count",
480 "diff btwn hi and lo for plane vs X cell per slice;slice delta plane;slicedelta X cell",
481 500,0.0,1000.0,250,0.0,500.0);
483 "diff btwn hi and lo for plane vs Y cell per slice;slice delta plane;slice delta Y cell",
484 500,0.0,1000.0,250,0.0,500.0);
486 "st. dev. of hit time per slice; slice sigma time [ns]; count",
487 500000,0.0,500000.0);
489 "Slice Duration (last hit time - first hit time);slice deltaT [ns];count",
493 "Plane vs. X cell for all hits (noise slice);noise slice plane;noise slice X cell",
494 1320,-20.0,1300.0,510,-10.0,500.0);
496 "Plane vs. Y cell for all hits (noise slice);noise slice plane;noise slice Y cell",
497 1320,-20.0,1300.0,510,-10.0,500.0);
499 "Time for all hits (noise slice); noise slice hit time [ns]; count",
500 12000,-50000.0,550000.0);
502 "avg hit plane vs. ave hit X cell (noise slice);avg noise slice plane;avg noise slice X cell",
503 660,-20.0,1300.0,255,-10.0,500.0);
505 "avg hit plane vs. ave hit Y cell (noise slice);avg noise slice plane;avg noise slice Y cell",
506 660,-20.0,1300.0,255,-10.0,500.0);
508 "avg hit time per slice (noise slice);avg noise slice time [ns]; count",
512 "Fraction of all hits in the noise slice; noise fraction; count",
515 fTres = tfs->
make<TH1F>(
"SliceHitTimeResolution",
516 "Computed timing resolution for hits in non-noise slices;non-noise hit time resolution [ns];count",
519 "Computed timing resolution for hits in the noise slice;noise hit time resolution [ns];count",
521 fPOT = tfs->
make<TH1F>(
"fPOT",
"total POT",1,0,1);
530 <<
"\n========== End Job ==========\n\n";
532 for(
int i = 1;
i <=
fNSvsVE->GetNbinsX(); ++
i) {
534 if(
fNSvsVE->GetBinContent(
i) != 0.0) {
584 for(
unsigned int i = 0;
i < slices->size(); ++
i)
587 std::vector< std::vector< cheat::NeutrinoEffPur > > ALLnuEffPur;
592 std::set< art::Ptr< simb::MCTruth > > NuListFLSHitALL;
593 std::set< art::Ptr< simb::MCTruth > > NuListCellHitALL;
594 std::set< art::Ptr< simb::MCTruth > > NuListCellHit;
596 std::set< art::Ptr< simb::MCTruth > > MCTListFLSHitALL;
597 std::set< art::Ptr< simb::MCTruth > > MCTListCellHitALL;
598 std::set< art::Ptr< simb::MCTruth > > MCTListCellHit;
614 for(
unsigned int i = 0;
i < slices->size(); ++
i) {
617 if( (*slices)[
i].IsNoise() ) {
625 std::vector< art::Ptr<caf::StandardRecord> > records = recordFMP.at(
i);
626 if (records.size() > 0 && recordFMP.isValid()) {
629 else { pass =
false; }
635 unsigned int IsNoise = 0;
655 double CXhi = -1.0e6;
657 double CYhi = -1.0e6;
660 float nu1X = 0.0, nu1Y = 0.0, nu1Z = 0.0, nu1E = -1.0;
661 int nu1pdg = 0, nu1CCNC = -1, nu1IntType = 0;
662 double nu1VisEsl = -1.0, nu1VisEtot = -1.0;
664 float nu2X = 0.0, nu2Y = 0.0, nu2Z = 0.0, nu2E = -1.0;
665 int nu2pdg = 0, nu2CCNC = -1, nu2IntType = 0;
666 double nu2VisEsl = -1.0, nu2VisEtot = -1.0;
668 float nu3X = 0.0, nu3Y = 0.0, nu3Z = 0.0, nu3E = -1.0;
669 int nu3pdg = 0, nu3CCNC = -1, nu3IntType = 0;
670 double nu3VisEsl = -1.0, nu3VisEtot = -1.0;
672 float nu4X = 0.0, nu4Y = 0.0, nu4Z = 0.0, nu4E = -1.0;
673 int nu4pdg = 0, nu4CCNC = -1, nu4IntType = 0;
674 double nu4VisEsl = -1.0, nu4VisEtot = -1.0;
676 float nu5X = 0.0, nu5Y = 0.0, nu5Z = 0.0, nu5E = -1.0;
677 int nu5pdg = 0, nu5CCNC = -1, nu5IntType = 0;
678 double nu5VisEsl = -1.0, nu5VisEtot = -1.0;
680 int pdgPartMost = 0, NumMatches = 0;
682 double nu1Eff = 0.0, nu1Pur = 0.0;
683 double nu2Eff = 0.0, nu2Pur = 0.0;
684 double nu3Eff = 0.0, nu3Pur = 0.0;
685 double nu4Eff = 0.0, nu4Pur = 0.0;
686 double nu5Eff = 0.0, nu5Pur = 0.0;
688 unsigned int IsMichelPrimary = 0;
689 unsigned int IsMichelPresent = 0;
690 unsigned int IsLateSlice = 0;
692 double PercentNoise = -1.0;
693 double PercentNoiseByPE = -1.0;
697 if((*slices)[
i].IsNoise()) IsNoise = 1;
698 NumHitsX = (*slices)[
i].NXCell();
699 NumHitsY = (*slices)[
i].NYCell();
703 if(!(*slices)[
i].IsNoise()) NSlice++;
710 double nx = 0.0, ny = 0.0;
713 for(
unsigned int j = 0;
j < (*slices)[
i].NCell(); ++
j) {
718 if((*slices)[
i].IsNoise()) {
734 double Tus = (h->
TNS())/1000.0;
745 if(h->
TNS() > Thi) Thi = h->
TNS();
746 if(h->
TNS() < Tlo) Tlo = h->
TNS();
752 if(h->
Cell() > CXhi) CXhi = h->
Cell();
753 if(h->
Cell() < CXlo) CXlo = h->
Cell();
759 if(h->
Cell() > CYhi) CYhi = h->
Cell();
760 if(h->
Cell() < CYlo) CYlo = h->
Cell();
769 if(nx == 0.0 && ny == 0.0) {
779 Psd =
sqrt(Psd/(nx+ny) - Pave*Pave);
780 Tsd =
sqrt(Tsd/(nx+ny) - Tave*Tave);
789 CXsd =
sqrt(CXsd/nx - CXave*CXave);
798 CYsd =
sqrt(CYsd/ny - CYave*CYave);
812 std::vector<cheat::NeutrinoEffPur> nuEffPur;
815 if(nuEffPur.size() > 0) {
818 if((*slices)[
i].IsNoise()) {
819 for(
unsigned int a = 0;
a < nuEffPur.size(); ++
a) {
835 for(
unsigned int a = 0;
a < nuEffPur.size(); ++
a) {
836 MCTListCellHitALL.insert(nuEffPur[
a].neutrinoInt);
837 if(nuEffPur[
a].neutrinoInt->NeutrinoSet())
838 NuListCellHitALL.insert(nuEffPur[
a].neutrinoInt);
842 if( !((*slices)[i].IsNoise()) ) {
843 for(
unsigned int a = 0;
a < nuEffPur.size(); ++
a) {
844 MCTListCellHit.insert(nuEffPur[
a].neutrinoInt);
845 if(nuEffPur[
a].neutrinoInt->NeutrinoSet())
846 NuListCellHit.insert(nuEffPur[
a].neutrinoInt);
854 if(nuEffPur[0].neutrinoInt->NeutrinoSet()) {
867 nu1Eff = nuEffPur[0].efficiency;
868 nu1Pur = nuEffPur[0].purity;
870 nu1VisEsl = nuEffPur[0].energySlice;
871 nu1VisEtot = nuEffPur[0].energyTotal;
873 if(nuEffPur.size() > 1) {
874 if(nuEffPur[1].neutrinoInt->NeutrinoSet()) {
876 nu = nuEffPur[1].neutrinoInt->GetNeutrino();
890 nu2Eff = nuEffPur[1].efficiency;
891 nu2Pur = nuEffPur[1].purity;
892 nu2VisEsl = nuEffPur[1].energySlice;
893 nu2VisEtot = nuEffPur[1].energyTotal;
895 if(nuEffPur.size() > 2) {
896 if(nuEffPur[2].neutrinoInt->NeutrinoSet()) {
898 nu = nuEffPur[2].neutrinoInt->GetNeutrino();
912 nu3Eff = nuEffPur[2].efficiency;
913 nu3Pur = nuEffPur[2].purity;
914 nu3VisEsl = nuEffPur[2].energySlice;
915 nu3VisEtot = nuEffPur[2].energyTotal;
917 if(nuEffPur.size() > 3) {
918 if(nuEffPur[3].neutrinoInt->NeutrinoSet()) {
920 nu = nuEffPur[3].neutrinoInt->GetNeutrino();
933 nu4Eff = nuEffPur[3].efficiency;
934 nu4Pur = nuEffPur[3].purity;
935 nu4VisEsl = nuEffPur[3].energySlice;
936 nu4VisEtot = nuEffPur[3].energyTotal;
938 if(nuEffPur.size() > 4) {
939 if(nuEffPur[4].neutrinoInt->NeutrinoSet()) {
941 nu = nuEffPur[4].neutrinoInt->GetNeutrino();
954 nu5Eff = nuEffPur[4].efficiency;
955 nu5Pur = nuEffPur[4].purity;
956 nu5VisEsl = nuEffPur[4].energySlice;
957 nu5VisEtot = nuEffPur[4].energyTotal;
961 NumMatches = nuEffPur.size();
964 const std::vector< const sim::Particle * > particles = BT->
HitsToParticle((*slices)[i].AllCells());
966 if(particles.size() > 0) {
967 pdgPartMost = particles[0]->PdgCode();
970 if(pdgPartMost == 11 && particles[0]->Process() ==
"Decay") IsMichelPrimary = 1;
972 for(
unsigned int a = 0;
a < particles.size(); ++
a) {
973 if(particles[
a]->
PdgCode() == 11 && particles[
a]->Process() ==
"Decay") IsMichelPresent = 1;
979 if(nuEffPur.size() > 0) {
980 if(nuEffPur[0].neutrinoInt->NParticles() > 0) {
983 if( (T-Tave) > 200.0) IsLateSlice = 1;
987 if(NumHitsX >= 3 && NumHitsY >= 3 &&
988 nu1Eff >= 0.9 && nu1Pur >= 0.9)
NgoodSL++;
994 double NNhitPE = 0.0;
996 for(
unsigned int hi = 0;
hi < (*slices)[
i].NCell(); ++
hi) {
1001 NhitPE += hit->
PE();
1005 NNhitPE += hit->
PE();
1010 PercentNoise = NNhit/Nhit;
1011 PercentNoiseByPE = NNhitPE/NhitPE;
1116 if((*slices)[
i].IsNoise()) {
1118 if(NumHitsX+NumHitsY == 0.0)
1119 std::cout <<
"\n\nEmpty Noise Slice! " << evt.
id() <<
"\n\n\n";
1125 NhitNoise += NumHitsX+NumHitsY;
1126 NhitTotal += NumHitsX+NumHitsY;
1132 fComp ->Fill(nu1Eff);
1133 fPur ->Fill(nu1Pur);
1143 NhitTotal += NumHitsX+NumHitsY;
1152 if(nu1Eff >= 0.9 && nu1Pur >= 0.9)
fNGvsVE->Fill(nu1VisEtot);
1172 for(
size_t p = 0;
p < parts->size(); ++
p) {
1173 for(
size_t t = 0;
t < truthAssn.at(
p).size(); ++
t) {
1175 MCTListFLSHitALL.insert(mcTru);
1177 NuListFLSHitALL.insert(mcTru);
1189 fNumNu->Fill(NuListFLSHitALL.size());
1190 fNumMCT->Fill(MCTListFLSHitALL.size());
1227 std::vector<cheat::NeutrinoEffPur>
results;
1230 for(
unsigned int i = 0;
i < nuEP.size(); ++
i) {
1231 if(nuEP[
i].
purity > 0.0) results.push_back(nuEP[
i]);
1236 if(results.size() > 0) {
1240 for(
unsigned int i = 1;
i < results.size(); ++
i) {
1243 results[
i] = results[
i-1];
1244 results[
i-1] =
temp;
double E(const int i=0) const
std::vector< int > nu2pdgVec
PDG code for second nu associated with slice.
std::vector< int > nu5pdgVec
PDG code for fifth nu associated with slice.
std::vector< double > nu4VisEtotVec
total visible energy for the fourth best matched nu in the entire event
std::vector< double > slEVec
slice.TotalGeV()
SubRunNumber_t subRun() const
std::string fCAFLabel
label for the process that made the standard records
back track the reconstruction to the simulation
std::vector< double > nu4VisEslVec
visible energy for the fourth best matched nu in the slice
double GetTimeRes(rb::CellHit const &cellhit)
std::vector< double > nu2VisEtotVec
total visible energy for the second best matched nu in the entire event
std::vector< unsigned int > IsMichelPresentVec
did a Michel e contribute to this slice?
int NumMCTCellHit
total number of MCTruths in non-noise slices counted by CellHits
std::vector< float > nu2YVec
y location for second nu associated with slice
std::vector< double > nu5VisEtotVec
total visible energy for the fifth best matched nu in the entire event
std::vector< int > pdgPartMostVec
the PDG code for the particle that contributed the most energy to this slice
std::vector< double > nu2VisEslVec
visible energy for the second best matched nu in the slice
void analyze(const art::Event &evt)
unsigned short Plane() const
std::vector< double > nu3EffVec
efficiency for third nu associated with slice
std::vector< double > PaveVec
average plane number for hits in a non-noise slice
std::vector< int > NumHitsYVec
number of Y hits per slice
const simb::MCParticle & Nu() const
std::vector< int > nu5IntTypeVec
Interaction type for fifth neutrino associated with slice.
std::vector< int > NumHitsXVec
number of X hits per slice
Vertical planes which measure X.
std::vector< float > nu4EVec
Energy for fourth nu associated with slice.
std::vector< double > CYhiVec
highest Y cell in the slice
std::vector< float > nu5YVec
y location for fifth nu associated with slice
std::vector< double > PercentNoiseVec
Percentage of slice that is true noise (by hit)
std::vector< unsigned int > IsMichelPrimaryVec
is this slice a Michel e? (the primary contribution by energy was from a michel e) ...
std::vector< double > PloVec
lowest plane in the slice
std::vector< int > nu4IntTypeVec
Interaction type for fourth neutrino associated with slice.
std::vector< double > TloVec
lowest plane in the slice
int NumMCTFLSHitALL
total number of MCTruths counted by FLSHits
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
recovalid::CAFCutter fCAFCutter
the CAFCutter helper class
DEFINE_ART_MODULE(TestTMapFile)
std::vector< float > nu3EVec
Energy for third nu associated with slice.
virtual void endSubRun(const art::SubRun &sr)
std::vector< float > nu1YVec
y location for primary nu associated with slice
std::vector< float > nu2ZVec
z location for second nu associated with slice
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
std::vector< double > nu3VisEtotVec
total visible energy for the third best matched nu in the entire event
std::vector< double > nu1EffVec
efficiency for primary nu associated with slice
Horizontal planes which measure Y.
std::vector< int > nu3pdgVec
PDG code for third nu associated with slice.
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
std::vector< double > nu5VisEslVec
visible energy for the fifth best matched nu in the slice
std::vector< double > PsdVec
st. dev. of plane number for hits in a non-noise slice
std::vector< float > nu2XVec
x location for second nu associated with slice
int NumMCTCellHitALL
total number of MCTruths counted by CellHits
std::vector< double > nu3VisEslVec
visible energy for the third best matched nu in the slice
std::vector< float > nu5XVec
x location for fifth nu associated with slice
std::vector< int > nu1CCNCVec
CCNC for primary nu associated with slice.
std::vector< double > nu1VisEslVec
visible energy for the best matched nu in the slice
std::vector< float > nu5ZVec
z location for fifth nu associated with slice
unsigned short Cell() const
std::vector< double > CXsdVec
st. dev. cell number for X hits in a non-noise slice
int InteractionType() const
bool passCuts(int cut, const caf::StandardRecord *sr)
TSpline3 hi("hi", xhi, yhi, 18,"0")
std::vector< double > nu1VisEtotVec
total visible energy for the best matched nu in the entire event
int NumNoise
number of noise slices
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
int NumNuFLSHitALL
number of nu counted by FLSHits
std::vector< double > nu5PurVec
purity for fifth nu associated with slice
std::vector< double > TsdVec
st. dev. time for hits in a non-noise slice
int fCutType
what cuts should we apply (see CAFCutter.h)
void push_back(Ptr< U > const &p)
int NumNuCellHit
number of nu in non-noise slices counted by CellHits
int NumSlice
number of non-noise slices
std::vector< double > PercentNoiseByPEVec
Percentage of slice that is true noise (weighted by PE)
std::vector< float > nu1EVec
Energy for primary nu associated with slice.
std::vector< double > nu2PurVec
purity for second nu associated with slice
std::vector< int > NumMatchesVec
number of MCTruth matches per slice
std::vector< float > nu5EVec
Energy for fifth nu associated with slice.
std::vector< float > nu2EVec
Energy for second nu associated with slice.
std::vector< double > CYaveVec
average cell number for Y hits in a non-noise slice
T get(std::string const &key) const
bool fApplyCAFCuts
should we apply CAF cuts?
double T(const int i=0) const
std::vector< double > nu3PurVec
purity for third nu associated with slice
std::vector< unsigned int > IsLateSliceVec
is the average time for this slice > 200 ns after the time of the neutrino interaction?
std::string fSlicerLabel
label for the process that made the slices
SlicerAna(fhicl::ParameterSet const &pset)
fvar< T > Phi(const fvar< T > &x)
EDAnalyzer(Table< Config > const &config)
std::vector< double > nu4PurVec
purity for fourth nu associated with slice
std::vector< double > nu2EffVec
efficiency for second nu associated with slice
std::vector< float > slPEVec
sum of PE for all hits in the slice
std::vector< int > nu5CCNCVec
CCNC for fifth nu associated with slice.
std::vector< double > CXaveVec
average cell number for X hits in a non-noise slice
bool fIsMC
is the file being processed MC?
std::vector< double > CXloVec
lowest X cell in the slice
EventNumber_t event() const
void reconfigure(const fhicl::ParameterSet &pset)
std::vector< int > nu3CCNCVec
CCNC for third nu associated with slice.
std::vector< int > nu1IntTypeVec
Interaction type for primary neutrino associated with slice.
std::vector< double > TaveVec
average time for hits in a non-noise slice
double Vx(const int i=0) const
std::vector< double > CYsdVec
st. dev. cell number for Y hits in a non-noise slice
A rawdata::RawDigit with channel information decoded.
T * make(ARGS...args) const
std::vector< int > nu2IntTypeVec
Interaction type for second neutrino associated with slice.
std::vector< int > nu3IntTypeVec
Interaction type for third neutrino associated with slice.
std::vector< unsigned int > IsNoiseVec
0 = non-noise slice, 1 = noise slice
std::vector< float > nu3YVec
y location for third nu associated with slice
int NumNuCellHitALL
number of nu counted by CellHits
std::vector< int > nu1pdgVec
PDG code for primary nu associated with slice.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
std::vector< float > nu3ZVec
z location for third nu associated with slice
std::vector< float > nu1ZVec
z location for primary nu associated with slice
std::vector< float > nu1XVec
x location for primary nu associated with slice
std::vector< cheat::NeutrinoEffPur > sortEffPur(const std::vector< cheat::NeutrinoEffPur > &nuEP)
std::vector< float > nu4YVec
y location for fourth nu associated with slice
double Vz(const int i=0) const
std::vector< float > nu4ZVec
z location for fourth nu associated with slice
std::vector< double > nu4EffVec
efficiency for fourth nu associated with slice
std::vector< double > CYloVec
lowest Y cell in the slice
std::vector< double > ThiVec
highest time in the slice
Class to help Slicer4D manage the collection of hits.
std::vector< float > nu4XVec
x location for fourth nu associated with slice
std::vector< double > PhiVec
highest plane in the slice
Helper class for Reco Validation modules.
std::vector< double > nu1PurVec
purity for primary nu associated with slice
std::vector< double > CXhiVec
highest X cell in the slice
std::vector< double > nu5EffVec
efficiency for fifth nu associated with slice
double totgoodpot
normalized by 10^12 POT
std::vector< float > nu3XVec
x location for third nu associated with slice
std::vector< int > nu4CCNCVec
CCNC for fourth nu associated with slice.
Event generator information.
std::vector< int > nu4pdgVec
PDG code for fourth nu associated with slice.
double Vy(const int i=0) const
float fTotalPOT
what is the POT in this file?
std::vector< int > nu2CCNCVec
CCNC for second nu associated with slice.
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.