20 #include "TMVA/Factory.h" 21 #include "TMVA/Tools.h" 22 #include "TMVA/Reader.h" 23 #include "TMVA/MethodKNN.h" 49 #include "Utilities/AssociationUtil.h" 78 void setTruthVars(
const std::vector<cheat::NeutrinoEffPur>& nuEP);
218 std::cerr <<
"\n\n\nWARNING: PID training mode for PDG code = " 220 <<
" is not currently supported. ABORTING...\n\n\n";
231 fdEdxVSbg[0][1] = (TH2F*)
infile.FindObjectAny(
"fdEdxVSbg_W0_X1");
232 fdEdxVSbg[1][0] = (TH2F*)
infile.FindObjectAny(
"fdEdxVSbg_W1_X0");
233 fdEdxVSbg[1][1] = (TH2F*)
infile.FindObjectAny(
"fdEdxVSbg_W1_X1");
259 fBPFTrackSignalTree = tfs->
make<TTree>(
"BPFTrackSignalTree",
"Tree for TMVA training with signal track level BPF info");
293 fBPFTrackBackgroundTree = tfs->
make<TTree>(
"BPFTrackBackgroundTree",
"Tree for TMVA training with background track level BPF info");
330 fBPFSliceTree = tfs->
make<TTree>(
"BPFSliceTree",
"Tree for TMVA training with slice level BPF info");
368 fPOT = tfs->
make<TH1D>(
"fPOT",
"Total POT",1,0,1);
444 float mostPure = 0.0;
445 unsigned int match = 0;
446 for(
unsigned int i = 0;
i < nuEP.size(); ++
i) {
447 if(nuEP[
i].
purity > mostPure && nuEP[
i].neutrinoInt->NeutrinoSet()) {
448 mostPure = nuEP[
i].purity;
453 if(mostPure > 0 && nuEP[match].neutrinoInt->NeutrinoSet()) {
503 for(
unsigned int i = 0;
i < slices->size(); ++
i) {
506 std::vector< std::vector< cheat::NeutrinoEffPur > > ALLnuEffPur;
512 for(
unsigned int islice = 0; islice < slices->size(); ++islice) {
515 if((*slices)[islice].IsNoise())
continue;
526 fMinX = (*slices)[islice].MinX();
527 fMinY = (*slices)[islice].MinY();
528 fMinZ = (*slices)[islice].MinZ();
529 fMaxX = (*slices)[islice].MaxX();
530 fMaxY = (*slices)[islice].MaxY();
531 fMaxZ = (*slices)[islice].MaxZ();
532 fNHitX = (*slices)[islice].NXCell();
533 fNHitY = (*slices)[islice].NYCell();
539 std::vector< art::Ptr< rb::Prong > > prongs;
540 if(prongFMP.isValid()) prongs = prongFMP.at(islice);
549 float longest = -1.0;
550 float bestPID = -1.0;
559 for(
unsigned int iprong = 0; iprong < prongs.size(); ++iprong) {
561 std::vector< art::Ptr< rb::Track > > tracks;
562 if(trackFMP.isValid()) tracks = trackFMP.at(iprong);
565 for(
unsigned int itrack = 0; itrack < tracks.size(); ++itrack) {
567 std::vector< art::Ptr< rb::FitSum > > fitsums;
568 if(fitsumFMP.isValid()) fitsums = fitsumFMP.at(itrack);
569 if(fitsums.size() != 1) abort();
571 float len = tracks[itrack]->TotalLength();
572 if(len > longest && fitsums[0]->PDG() ==
fPIDmode) {
588 for(
unsigned int iprong = 0; iprong < prongs.size(); ++iprong) {
591 std::vector< art::Ptr< rb::Track > > tracks;
592 if(trackFMP.isValid()) tracks = trackFMP.at(iprong);
598 float chi2Tmu = 1.0e9;
599 float chi2Tpi = 1.0e9;
600 float chi2Tpr = 1.0e9;
603 for(
unsigned int itrack = 0; itrack < tracks.size(); ++itrack) {
606 std::vector< art::Ptr< rb::FitSum > > fitsums;
607 if(fitsumFMP.isValid()) fitsums = fitsumFMP.at(itrack);
612 if(fitsums.size() != 1) abort();
622 tempH = fitsums[0]->Chi2(0)+fitsums[0]->Chi2(1);
623 tempA = fitsums[0]->Chi2(2)+fitsums[0]->Chi2(3);
624 tempN = fitsums[0]->Ndof(0)+fitsums[0]->Ndof(1)-4;
627 if(fitsums[0]->PDG() ==
fPIDmode && tempN > 0.0) {
631 fLength = tracks[itrack]->TotalLength();
635 if(fitsums[0]->PDG() == 13 && tempN > 0.0) chi2Tmu = tempH/tempN+tempA/tempN;
636 if(fitsums[0]->PDG() == 211 && tempN > 0.0) chi2Tpi = tempH/tempN+tempA/tempN;
637 if(fitsums[0]->PDG() == 2212 && tempN > 0.0) chi2Tpr = tempH/tempN+tempA/tempN;
645 std::vector<double>
dE;
646 std::vector<double>
dx;
647 std::vector<double>
W;
648 std::vector<double> BG;
649 std::vector<double>
s;
656 for(
unsigned int a = 0;
a < dE.size(); ++
a) {
661 if(W[
a] > 0.0) w = 1;
662 if(dx[
a] > 8.0) x = 1;
669 if(content < 1.0
e-9) content = 1.0e-9;
670 LLtemp +=
log(content);
694 prongPick = (
int)iprong;
695 fPmu = fitsums[0]->FourMom().P();
696 fDirZmu = tracks[itrack]->Dir().Z();
711 if(tracks.size() == 3) {
719 const std::vector< const sim::Particle * > particles = BT->
HitsToParticle(prongs[iprong]->AllCells());
754 for(
unsigned int p = 0;
p < prongs.size(); ++
p) {
757 if(prongs[
p]->Is2D())
continue;
764 for(
unsigned int h = 0;
h < prongs[
p]->NCell(); ++
h) {
770 for(
unsigned int hpr = 0; hpr < hits3Dprongs.
NCell(); ++hpr) {
772 if((*pronghit) == (*pronghit3D)) {
779 if(!onlist) hits3Dprongs.
Add(prongs[
p]->Cell(
h));
784 if((
int)
p == prongPick)
continue;
799 for(
unsigned int h = 0;
h < prongs[
p]->NCell(); ++
h) {
802 float W = prongs[
p]->W(chit);
806 if(rhit.IsCalibrated()) E3Dpr += rhit.GeV();
825 for(
unsigned int hsl = 0; hsl < (*slices)[islice].NCell(); ++hsl) {
831 for(
unsigned int hpr = 0; hpr < hits3Dprongs.
NCell(); ++hpr) {
833 if((*pronghit) == (*slicehit)) {
841 hitsRemain.
Add((*slices)[islice].Cell(hsl));
842 sumPE += slicehit->
PE();
858 float muonTrackHadE = 0.0;
863 std::vector< art::Ptr< rb::Track > > tracks;
864 if(trackFMP.isValid()) tracks = trackFMP.at(prongPick);
868 for(
unsigned int itrack = 0; itrack < tracks.size(); ++itrack) {
870 std::vector< art::Ptr< rb::FitSum > > fitsums;
871 if(fitsumFMP.isValid()) fitsums = fitsumFMP.at(itrack);
872 if(fitsums.size() != 1) abort();
875 if(fitsums[0]->PDG() != 13)
continue;
898 if(KtrackFMP.isValid()) {
900 std::vector< art::Ptr<rb::Track> > KalTracks = KtrackFMP.at(islice);
902 float bestRemid = -5.0;
904 if(!remidFMP.isValid())
continue;
906 for(
unsigned int ikt = 0; ikt < KalTracks.size(); ++ikt) {
907 std::vector< art::Ptr<remid::ReMId> > remids = remidFMP.at(ikt);
909 if(remids.size() != 1)
continue;
910 if(remids[0]->Value() > bestRemid) bestRemid = remids[0]->Value();
922 const std::vector< const sim::Particle * > particles = BT->
HitsToParticle(prongs[prongPick]->AllCells());
924 if(particles.size() > 0) pdg =
abs(particles[0]->
PdgCode());
double E(const int i=0) const
back track the reconstruction to the simulation
unsigned int fTrainMode
training mode (energy estimator or muon PID?)
std::string fXMLFile
location of the XML file for the muon PID
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
float fLength
track length
std::string fHistoFile
location of dE/dx log-likelihood histogram
double fTotalPOT
total POT counted for all events in this job
TTree * fBPFSliceTree
TTree with slice level info (for energy estimator training)
murem::TrackCleanUpAlg fTrackCleanUpAlg
object to determine hadronic energy on the muon track
float fEsq
true energy squared (used to estiamte st. dev.)
float fDChi2Tmupi
delta chi2 total muon - pion
const simb::MCParticle & Nu() const
int fBestBPFmuPIDpdg
pdg code for the particle that contributed the most energy to the track with the best BPF muon PID ...
int fnuPDG
PDG code for the best matched nu (from truth)
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
void analyze(const art::Event &evt)
A collection of associated CellHits.
int fCCNC
is this CC or NC (from truth)
float fMaxY
maximum hit Y value in the slice
void computeDEDX(art::Ptr< rb::Track > &track, int pdg=13)
Compute dE/dx for all cells along this track and the fitsum object that went with it...
float fSumPE
sum of the PE for all hits not on the muon track
int fNHitY
number of Y-view hits in the slice
DEFINE_ART_MODULE(TestTMapFile)
float fDChi2Tpipr
delta chi2 total pion - proton
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...
bpfit::dEdxCalculator fdEdxCalc
helper class for computing dEdx values
int fPIDmode
make the training TTree using the muon, pion, or proton fit
float fMaxX
maximum hit X value in the slice
float fLengthRatio
ratio of track to longest track in the slice
float fN3Dprongs
total number of 3D prongs
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
TH1D * fPOT
Total POT accumulated for all events processed in this job.
const art::ProductToken< sumdata::POTSum > fPOTSumToken
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Calculates dE/dx in cells passed through by a track.
float fChi2H
chi^2 hits value per DOF
std::string fProngInstance
instance label for the prongs to be used
unsigned short Cell() const
float fMinZ
minimum hit Z value in the slice
TTree * fBPFTrackBackgroundTree
TTree with track level background info (for PID training)
int InteractionType() const
TMVA::Reader fMuonPID
TMVA reader object for the muon PID.
int fNHitX
number of X-view hits in the slice
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
void push_back(Ptr< U > const &p)
float fE3Dprongs
sum of hit energies in 3D prongs (excluding the "muon" prong) - this does NOT account for dead materi...
T get(std::string const &key) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Near Detector in the NuMI cavern.
float fDirZmu
muon track dirZ from the track selected to be the muon
float fEremain
sum of remaining slice energy not in 3D prongs PLUS extra energy added back in from TrackCleanUpAlg ...
float fTrackPDG
PDG code for the primary track particle.
double fdxTOL
lower limit for dx values to be used in dE/dx calculation
EDAnalyzer(Table< Config > const &config)
Perform a "2 point" Hough transform on a collection of hits.
virtual void endSubRun(const art::SubRun &sr)
TH2F * fdEdxVSbg[2][2]
histo of computed dE/dx for cell hits (vs. beta*gamma)
void setTruthVars(const std::vector< cheat::NeutrinoEffPur > &nuEP)
std::string fProngLabel
label for module that made the prongs
std::string fTrackLabel
label for module that made the tracks
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
float fdEdxLL
dEdx LL value
float fBestBPFmuPID
best BPF muon PID value
float fPmu
reconstructed muon momentum from the track selected to be the muon
A rawdata::RawDigit with channel information decoded.
T * make(ARGS...args) const
float fHitRatio
ratio of track hits to prong hits (used as a track/shower discriminator)
float fMaxZ
maximum hit Z value in the slice
float fMinX
minimum hit X value in the slice
float fEhadMu
Hadronic energy contamination on the muon track as computed by TrackCleanUpAlg.
Service to store calibration data products (CDP) in the SQLite3 metadatabase of a file...
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
BPFTmvaTrainer(fhicl::ParameterSet const &pset)
const art::ProductToken< std::vector< rb::Cluster > > fSlicerToken
float fDChi2Tmupr
delta chi2 total muon - proton
int fIntType
neutrino interaction type (from truth)
void reconfigure(const fhicl::ParameterSet &pset)
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
double totgoodpot
normalized by 10^12 POT
Event generator information.
TTree * fBPFTrackSignalTree
TTree with track level signal info (for PID training)
float fBestRemid
best remid value
float ExtraEOnTrackInGeV(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
ProductToken< T > consumes(InputTag const &)
float fChi2A
chi^2 angle value per DOF
void getResults(std::vector< double > &dEtemp, std::vector< double > &dxtemp, std::vector< double > &Wtemp, std::vector< double > &bgtemp, std::vector< double > &stemp, double dxTol)
Return the results of the dE/dx calculation.
Encapsulate the geometry of one entire detector (near, far, ndos)
float fMinY
minimum hit Y value in the 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.