17 #include "Utilities/AssociationUtil.h" 111 nHits[0] = nHits[1] = 0;
112 nOnlyHits[0] = nOnlyHits[1] = 0;
143 TreeS = tfs->
make<TTree>(
"fTreeS",
"kNN Input Variables for Signal");
153 TreeB = tfs->
make<TTree>(
"fTreeB",
"kNN Input Variables for Background");
168 DedxSample = tfs->
make<TTree>(
"fDedxSample",
"Variables for Creating Dedx Sample Histograms");
181 ScatSample = tfs->
make<TTree>(
"fScatSample",
"Variables for Creating Scattering Sample Histograms");
211 for(
unsigned int i = 0;
i<slicecol->size();++
i){
214 for(
unsigned int ihit = 0; ihit<slice->
NCell();++ihit){
229 while(!((dibmask>>iD)&1)){
245 unsigned int mostEffSlice = 0;
248 for(
unsigned int iSlice = 0;iSlice<slicelist.
size();++iSlice){
249 if(slicelist[iSlice]->IsNoise()){
continue; }
252 if(nuint.size() > 0 && nuint[0].efficiency > bestEff){
253 bestEff = nuint[0].efficiency;
254 mostEffSlice = iSlice;
258 fvisE = nuint[0].energySlice;
262 else if(nu.
Mode() == 0){nuintType = 0;}
263 else if(nu.
Mode() == 1){nuintType = 1;}
264 else if(nu.
Mode() == 2){nuintType = 2;}
265 else if(nu.
Mode() == 3){nuintType = 3;}
271 for(
unsigned int iSlice = 0; iSlice<slicelist.
size(); ++iSlice){
273 if(iSlice != mostEffSlice){
continue; }
275 if(bestEff == 0){
continue; }
277 if(slicelist[iSlice]->IsNoise()){
continue; }
279 if(nuintType == 6){
continue; }
282 std::vector<art::Ptr<rb::Track> > tracks = fndmnytrk.at(iSlice);
284 if(tracks.size() == 0){
continue; }
286 frecoE = slicelist[iSlice]->TotalGeV();
289 std::sort(tracks.begin(),tracks.end(),
SortByHits);
295 std::set<int> allowedDaughters;
296 allowedDaughters.insert(-11);
297 allowedDaughters.insert(11);
304 std::set<int> primTrkIDs;
306 for(std::set<int>::iterator
it = trackIDs.begin();
it != trackIDs.end();++
it){
311 std::vector<cheat::ParticleEffPur> pep = bt->
TracksToParticles(tracks, primTrkIDs, allowedDaughters, slicelist[iSlice]->AllCells());
317 for(
unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
319 fpdg = pep[itrk].pdg;
323 if(
fpdg == 0){
continue; }
331 fnPlane = tracks[itrk]->ExtentPlane();
345 if(tracks[itrk]->TotalLength() > 6500){
continue; }
348 if(!tracks[itrk]->Is3D()){
continue; }
354 for(
unsigned int idedx = 0; idedx<pid->
NDedx(); ++idedx){
366 for(
unsigned int iscat = 0; iscat<pid->
NScatters();++iscat){
370 TVector3 trkstart = tracks[itrk]->Start();
373 for(
unsigned int itraj = 0; itraj<tracks[itrk]->NTrajectoryPoints()-1; ++itraj){
374 TVector3
p0 = tracks[itrk]->TrajectoryPoint(itraj);
375 TVector3
p1 = tracks[itrk]->TrajectoryPoint(itraj+1);
376 lastdist = (p1-
p0).
Mag();
399 for(
unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
400 if(
abs(pep[itrk].
pdg) == 13){ bestMu = itrk; }
426 double longestPi = -1;
427 int longestTrkIdx = -1;
428 double longestTrk = -1;
431 for(
unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
432 double trklen = tracks[itrk]->TotalLength();
433 if(
abs(pep[itrk].
pdg) == 211 && trklen > longestPi){
437 if(trklen > longestTrk){
438 longestTrkIdx = itrk;
455 fpdg = pep[bestPi].pdg;
463 else if(longestTrkIdx >= 0){
474 fpdg = pep[longestTrkIdx].pdg;
510 bool contained =
false;
516 unsigned int slcncellsfromedge = 999999999;
517 for(
unsigned int ihit = 0; ihit < slice->
NCell(); ++ihit){
521 if(cellsfromedge < slcncellsfromedge){
522 slcncellsfromedge = cellsfromedge;
527 unsigned int firstpl = trk->
MinPlane();
528 unsigned int lastpl = trk->
MaxPlane();
530 TVector3 dirbwk = -trk->
Dir();
535 if(slcncellsfromedge > 1){
536 unsigned int npltofront = firstpl - (
fdibfirst-1)*64;
537 unsigned int npltoback = (
fdiblast*64) - 1 - lastpl;
539 if(npltofront > 1 && npltoback > 1){
542 if(fwdcell > 10 && bwkcell > 10){
551 if(slcncellsfromedge > 1){
557 if(firstpl > 1 && lastpl < 212){
562 if(fwdcell > 4 && bwkcell > 8 && trk->
Start().Z() < 1150){
565 if(trk->
Stop().Z() < 1275){
572 double mindist = 99999999;
573 unsigned int mintrajid = 0;
576 if(dist > 0 && dist < mindist){
582 TVector3 trandir = trk->
Dir();
583 if(mintrajid < trk->NTrajectoryPoints() -2){
double E(const int i=0) const
size_t NTrajectoryPoints() const
back track the reconstruction to the simulation
double Dedx(unsigned int i) const
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
unsigned int NScatters() const
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
bool fCreateTrainTrees
Create kNN training tree?
bool fCreateSampleTrees
Create trees for sample histograms?
const sim::Particle * TrackIDToMotherParticle(int const &id) const
std::vector< NeutrinoEffPur > SliceToNeutrinoInteractions(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of neutrino interactions...
unsigned short Plane() const
double ScatSeparation() const
Return the scattering separation variable used as an input to the kNN that determines the pid value...
const simb::MCParticle & Nu() const
art::ServiceHandle< geo::Geometry > fgeom
unsigned int Ncells() const
Number of cells in this plane.
std::string fTrackLabel
Where to find recondtructed tracks.
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
int ProjectedLiveCellsToEdge(TVector3 vertex, TVector3 dir)
recommended projection to use
bool DedxVertex(unsigned int i) const
unsigned short Cell() const
Track finder for cosmic rays.
Far Detector at Ash River, MN.
TTree * TreeS
Tree for holding signal input variables for kNN training.
void push_back(Ptr< U > const &p)
double YPositionAtTransition(TVector3 vertex, TVector3 dir)
int NMeasurementPlanes() const
T get(std::string const &key) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
virtual TVector3 Dir() const
Unit vector describing prong direction.
TVector3 ScatLocation(unsigned int i) const
int FullNDProjectedCells(TVector3 vertex, TVector3 dir)
Near Detector in the NuMI cavern.
double CosScat(unsigned int i) const
Collect Geo headers and supply basic geometry functions.
TTree * ScatSample
Tree for holding scattering variables needed for making signal sample histograms. ...
int NVertexPlanes() const
bool DedxUsed(unsigned int i) const
EDAnalyzer(Table< Config > const &config)
ReMIdTrain(fhicl::ParameterSet const &pset)
TTree * TreeB
Tree for holding background input variables for kNN training.
unsigned int NDedx() const
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
art::ServiceHandle< geo::LiveGeometry > flivegeom
std::string fSliceLabel
Where to find reconstructed slices.
T * make(ARGS...args) const
double DedxLocation(unsigned int i) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
TTree * DedxSample
Tree for holding dE/dx variables needed for making signal sample histograms.
std::string fPidLabel
Where to find the pids.
void analyze(art::Event const &evt)
const std::set< int > GetTrackIDList() const
Get all G4 track ids present in the event.
bool ContainedEvent(art::Ptr< rb::Cluster > slice, art::Ptr< rb::Track > trk)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
art::ServiceHandle< nova::dbi::RunHistoryService > frh
T min(const caf::Proxy< T > &a, T b)
TVector3 Stop() const
Position of the final trajectory point.
int GoodDiBlockMask(int subrun=-1, bool reload=false)
Event generator information.
std::vector< ParticleEffPur > TracksToParticles(const std::vector< const rb::Track * > &tracks, const std::set< int > &trkIDs, const std::set< int > &allowedDaughters, const std::vector< const rb::CellHit * > &allHits) const
Given a collection of reconstructed tracks, this method finds the best match for each reco track to t...
double DedxSeparation() const
Return the dE/dx separation variable used as an input to the kNN that determines the pid value...
Encapsulate the geometry of one entire detector (near, far, ndos)
static bool SortByHits(const art::Ptr< rb::Track > &a, const art::Ptr< rb::Track > &b)