21 #include "NovaDAQConventions/DAQConventions.h" 28 #include "Utilities/AssociationUtil.h" 33 #include "TMVA/Reader.h" 118 produces< std::vector<cosrej::CosRejObj> >();
119 produces< art::Assns<cosrej::CosRejObj, rb::Cluster> >();
120 produces< std::vector<cosrej::TrkCntObj> >();
121 produces< art::Assns<cosrej::TrkCntObj, rb::Track> >();
239 fReaderUCTuned->BookMVA(
"UCTuned_BDTA",pidpath+
"UCTuned_prod3_BDTA_700.weights.xml");
249 fReaderTAP1->BookMVA(
"TAP1_BDTA",pidpath+
"TAP1_BDTA.weights.xml");
259 fReaderTAP2->BookMVA(
"TAP2_BDTA",pidpath+
"TAP2_BDTA.weights.xml");
269 fReaderTAP3->BookMVA(
"TAP3_BDTA",pidpath+
"TAP3_BDTA.weights.xml");
286 fReaderProd5_Per1 -> BookMVA(
"BDT", pidpath+
"TMVA_MProd5_Kal_FHC_Per1_BDT.weights.xml");
303 fReaderProd5_Per2 -> BookMVA(
"BDT", pidpath+
"TMVA_MProd5_Kal_FHC_Per2_BDT.weights.xml");
320 fReaderProd5_FHC -> BookMVA(
"BDT", pidpath+
"TMVA_MProd5_Kal_FHC_High_BDT.weights.xml");
337 fReaderProd5_RHC -> BookMVA(
"BDT", pidpath+
"TMVA_MProd5_Kal_RHC_High_BDT.weights.xml");
349 if(slicevec->empty()) {
356 std::unique_ptr< std::vector<cosrej::CosRejObj> > outputObjects(
new std::vector<cosrej::CosRejObj>);
358 std::unique_ptr< std::vector<cosrej::TrkCntObj> > outTrkObjects(
new std::vector<cosrej::TrkCntObj>);
370 mf::LogError(
"CosRej") <<
"Spill Data not found, aborting without horn current information";
373 isRHC = spillPot->
isRHC;
379 for(
size_t i = 0;
i < slicevec->size(); ++
i){
389 unsigned int bestIdx = 999;
392 double fwddist = 0, bakdist = 0, scatt = 0;
393 int fwdcell = 0, bakcell = 0;
394 bool badtrack =
false;
396 const std::vector< art::Ptr<rb::Track> > cosTracks = costrackAssnList.at(
i);
401 if(!costrackAssnList.isValid()) {
402 mf::LogError(
"CosRej") <<
"CosRej: No Cosmic Tracks! Do not go any further!";
403 outputObjects->push_back(cosrejobj);
408 if(cosTracks.size()==1) {
411 if(angleOfBestTrack>-6&&angleOfBestTrack<2) cosrejobj.
SetAngleCos(angleOfBestTrack);
415 TVector3 dirbak = -track->
Dir().Unit();
436 double slope = 0, chisq = 0, chisqdiff = 0;
454 const TVector3 trueDir = parts[0]->Momentum().Vect().Unit();
460 if(!trackAssnList.isValid()) {
462 outputObjects->push_back(cosrejobj);
467 const std::vector< art::Ptr<rb::Track> > sliceTracks = trackAssnList.at(
i);
472 for(
size_t iTrack = 0; iTrack < sliceTracks.size(); ++iTrack)
479 outTrkObjects->push_back(trkcntobj);
485 double tdx = track->
Dir().X();
486 double tdy = track->
Dir().Y();
487 double tdz = track->
Dir().Z();
490 if(!(tdx>=-1&&tdx<=1&&tdy>=-1&&tdy<=1&&tdz>=-1&&tdz<=1)){
491 outTrkObjects->push_back(trkcntobj);
498 TVector3 dirbak = -track->
Dir().Unit();
535 outTrkObjects->push_back(trkcntobj);
540 outputObjects->push_back(cosrejobj);
549 if(sliceTracks.size() == 0) {
550 outputObjects->push_back(cosrejobj);
555 outputObjects->push_back(cosrejobj);
560 else if(!remidVec.isValid()) {
562 outputObjects->push_back(cosrejobj);
568 for(
size_t j = 0;
j < sliceTracks.size(); ++
j){
571 if(!track->
Is3D() || trackRemid->
Value() < 0)
continue;
574 if(trackRemid->
Value() < 0.5) nhadtracks++;
579 if(cosrejobj.
NTracks3D()>0 && bestIdx!=999) {
587 unsigned int bestIdx2nd = 999;
588 double highestPID2nd = -999.0;
589 double bestTrkLen2nd = 0.0;
590 for(
size_t j = 0;
j < sliceTracks.size(); ++
j){
591 if(
j == bestIdx)
continue;
593 if(trackRemid2nd->
Value() < highestPID2nd)
continue;
594 if(trackRemid2nd->
Value() == highestPID2nd){
595 if(sliceTracks[
j]->TotalLength() < bestTrkLen2nd)
continue;
597 highestPID2nd = trackRemid2nd->
Value();
598 bestTrkLen2nd = sliceTracks[
j]->TotalLength();
605 if((angleOf2ndBestTrack>-6) && (angleOf2ndBestTrack<2)) cosrejobj.
SetAngleKal2nd(angleOf2ndBestTrack);
606 TVector3 firstDir = track->
Dir().Unit();
607 TVector3 secondDir = track2nd->
Dir().Unit();
608 double angleBtwKalTrks = (firstDir*secondDir);
623 const TVector3 trueDir = parts[0]->Momentum().Vect().Unit();
629 double tdx = track->
Dir().X();
630 double tdy = track->
Dir().Y();
631 double tdz = track->
Dir().Z();
633 if(!(tdx>=-1&&tdx<=1&&tdy>=-1&&tdy<=1&&tdz>=-1&&tdz<=1)) badtrack =
true;
635 if(badtrack)
continue;
639 TVector3 dirbak = -track->
Dir().Unit();
679 double slope = 0, chisq = 0, chisqdiff = 0;
693 if(angleOfBestTrack>-6&&angleOfBestTrack<2) cosrejobj.
SetAngleKal(angleOfBestTrack);
700 float fslcCalEPerNHit = 0.0;
701 float fhadEPerNHit = 0.0;
702 int frun = evt.
run();
705 if(fmCVN.isValid()) {
706 const std::vector<art::Ptr<cvn::Result>>
results = fmCVN.at(
i);
707 if(!results.empty()) {
720 TMVAvarsUCTuned[7] = fhadEPerNHit;
721 TMVAvarsUCTuned[8] = fslcCalEPerNHit;
724 TMVAvarsUCTuned[11] = cosrejobj.
Scatt();
737 if (frun > 10000 && frun < 17300) {
739 }
else if(frun > 17300 && frun < 20753) {
741 }
else if (frun > 20752) {
746 mf::LogError(
"CosRej") <<
"CVN Data Product " <<
fCVNLabel <<
" not found, needed for CosRej2019, aborting.";
756 TMVAvarsProd5[4] =
cos( track->
Dir().Y() );
758 TMVAvarsProd5[6] =
sqrt( 1 - (TMVAvarsProd5[1]*TMVAvarsProd5[1]) );
763 if ( frun <= 17139 ) {
765 }
else if ( frun <= 19746 ) {
791 outputObjects->push_back(cosrejobj);
796 evt.
put(std::move(outputObjects));
797 evt.
put(std::move(assoc));
798 evt.
put(std::move(outTrkObjects));
799 evt.
put(std::move(trkAssoc));
void SetKalChisq(float kalchisq)
chisq result associated with slope fit. Mostly not useful
T max(const caf::Proxy< T > &a, T b)
void SetCDirScore(float cdirscore)
chisq diff between assuming forwards going and backwards going for Cosmic Track
size_t NTrajectoryPoints() const
bool isRHC
is the beam in antineutrino mode, aka RHC
void SetKalFwdDist(float kalfwddist)
distance (dist,cm) of best track from end to det wall projected forwards (Fwd), based on Kalman Track...
back track the reconstruction to the simulation
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
void SetKFitSpeed(float kfitspeed)
inverse speed (ns/cm) of track fit by TimingFit module (negative indicates backwards track) for best ...
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
void SetTrkLenInCat(double trklenincat)
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
fvar< T > fabs(const fvar< T > &x)
void SetKalSlope(float kalslope)
slope of timing fit; a positive slope indicates track is properly directed
double ProjectedSteelDist(TVector3 vertex, TVector3 dir)
TMVA::Reader * fReaderProd5_FHC
void SetTrackE2nd(float tracke2nd)
Calorimetric energy of 2nd highest ReMId kalman track.
void SetKalBakSteel(float kalbaksteel)
distance (dist,cm) through muon catcher traveled, when projected forwards to detector edge...
double DistToEdgeXY(TVector3 vertex)
void SetCosBakAir(float cosbakair)
distance (dist,cm) through air traveled, when projected backwards to detector edge. For ND only, muon catcher air gap. NYI - just use ytrans
void SetKalChiDiff(float kalchidiff)
chisq difference of assuming + or - c as slope. A positive chisqdiff means the track is backwards...
void SetOldCosPID(float oldcospid)
PID for contained cosmic rejection (Old Alg, Prod 3 training)
void SetKScoreDiff(float kscorediff)
chisq diff between best directed fit and free fit (may indicate if things are wonky) for best ReMId K...
void SetCosFwdCell(int cosfwdcell)
cells in that line segment instead of distance of line segment
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
CosRej(fhicl::ParameterSet const &pset)
Store results of TimingFitAlg.
void SetAngleKal(float angleb)
cos of angle of best Kalman track wrt beam
void SetTrkBakCellND(int trkbakcellnd)
void SetNHadTracks(int nhadtracks)
number of 3D Kalman tracks with ReMId value of < 0.4 (not mu-like, or hadronic-like) ...
void SetCosFwdDist(float cosfwddist)
distance (dist,cm) of track from end to det wall projected forwards (Fwd), based on Cosmic Tracker (C...
void SetTrackLenKal2nd(float tracklenkal2nd)
Total length of 2nd highest ReMId kalman track.
void SetScatt(float scatt)
basic scattering variable; sum of all angular deviations between best kal track trajectory points / t...
void SetCosFwdAir(float cosfwdair)
distance (dist,cm) through air traveled, when projected forwards to detector edge. For ND only, muon catcher air gap. NYI
void SetVtxDist(float vtxdist)
Algorithm to determine the direction of tracks using timing.
void SetCosYPosAtTrans(float cosyposattrans)
Y position at transition to muon catcher, for determining if track went through air gap (ND only) ...
void SetTrkBakAir(float trkbakair)
distance (dist,cm) through muon catcher traveled, when projected forwards to detector edge...
std::string fNuMIBeamLabel
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
void SetTrkFwdAir(float trkfwdair)
cosmic track distance projected backwards from track start
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
void SetMinDist(float mindist)
minimum projected distance, forwards or backwards, of ANY Kalman track with more than 15 hits ...
void SetTrkLenInAct(double trkleninact)
void SetKalBakCellND(int kalbakcellnd)
note that fwd<->bak values are interchanged if track direction is backwards, including MC...
void SetCosBakCell(int cosbakcell)
note that fwd<->bak values are interchanged if track direction is backwards
void SetUnconTunedCosPID(float uncontunedcospid)
Tuned PID for uncontained cosmic rejection (TMVA) - Jose's.
TimingFitResult HoughFit(const rb::Track *track) const
int ProjectedLiveCellsToEdge(TVector3 vertex, TVector3 dir)
recommended projection to use
void SetNTracks3D(int ntracks3d)
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
void SetTrkYPosAtTrans(float trkyposattrans)
void SetStartAct(float startact)
measure of activity near best Kalman track start position; start pos is a good measure of vertex ...
double MCFrontZ() const
Z position at which the muon catcher starts.
void produce(art::Event &evt)
ProductID put(std::unique_ptr< PROD > &&product)
virtual double TotalLength() const
Length (cm) of all the track segments.
float TrkEPerNHit() const
bool FScatterEstim(art::Ptr< rb::Track > track, bool &fromTrkEnd, double &angleExt, double &angleSigma, double &angleMax, double &angleSum)
double ProjectedLiveDistanceToEdge(TVector3 vertex, TVector3 dir)
void SetKalFwdCell(int kalfwdcell)
cells in that line segment instead of distance of line segment
Far Detector at Ash River, MN.
void SetAngleBtwKalTrks(float anglebtw)
cos of angle of 2nd best Kalman track wrt beam
void SetTrkFwdCellND(int trkfwdcellnd)
void SetCosChiDiff(float coschidiff)
same for Cosmic track
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
void SetTrkFwdSteel(float trkfwdsteel)
distance (dist,cm) through air traveled, when projected forwards to detector edge. For ND only, muon catcher air gap. NYI
double DistToEdgeZ(TVector3 vertex)
void SetCosThetaTrue(float cosqtrue)
cosine of angle between true particle (bt most contributing) and reco track dir; use to tell if track...
void getFits(const art::Ptr< rb::Track > track, double &slope, double &chisq, double &chisqdiff)
void SetEndDist(float enddist)
double YPositionAtTransition(TVector3 vertex, TVector3 dir)
float getScatt(const art::Ptr< rb::Track > track)
void SetCosBakCellND(int cosbakcellnd)
note that fwd<->bak values are interchanged if track direction is backwards, including MC...
void SetTrkFwdCell(int trkfwdcell)
float getAngle(TVector3 trackdir)
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
virtual TVector3 Dir() const
Unit vector describing prong direction.
int FullNDProjectedCells(TVector3 vertex, TVector3 dir)
Near Detector in the NuMI cavern.
TMVA::Reader * fReaderTAP3
void SetCFitSpeed(float cfitspeed)
inverse speed (ns/cm) of track fit by TimingFit module (negative indicates backwards track) for Cosmi...
art::ServiceHandle< geo::LiveGeometry > livegeom
void SetTrackLenKal(float tracklenkal)
Total length of kalman track.
void SetKalYPosAtTrans(float kalyposattrans)
Y position at transition to muon catcher, for determining if track went through air gap (ND only) ...
Perform a "2 point" Hough transform on a collection of hits.
void SetCosFwdSteel(float cosfwdsteel)
distance (dist,cm) through muon catcher traveled, when projected forwards to detector edge...
void SetMinCell(int mincell)
minimum projected number of cell, forwards or backwards, of ANY Kalman track with more than 15 hits ...
TMVA::Reader * fReaderProd5_RHC
std::string fCosmicTrackLabel
void SetCScoreDiff(float cscorediff)
chisq diff between best directed fit and free fit (may indicate if things are wonky) for Cosmic Track...
virtual double DistanceFromStart(double z) const
void SetKDirScore(float kdirscore)
chisq diff between assuming forwards going and backwards going for best ReMId Kalman track ...
void SetKalThetaTrue(float kalqtrue)
cosine of angle between true particle (bt most contributing) and reco track dir; use to tell if track...
void SetPDGBest(int pdgbest)
pdg code of true most contributing particle to reco track hits. Calculate for best Kalman track...
TMVA::Reader * fReaderProd5_Per2
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
void SetFScattSum(float fscattsum)
sum of scattering angles (Fernanda's)
std::vector< TrackIDE > HitsToTrackIDE(const std::vector< const rb::CellHit * > &hits, bool useBirksE=false) const
Returns vector of TrackIDE structs contributing to the given collection of hits.
Cosmic Rejection PIDs for Numu analysis.
float TrackLenKal() const
float TMVAvarsUCTuned[74]
void SetKalBakAir(float kalbakair)
distance (dist,cm) through air traveled, when projected forwards to detector edge. For ND only, muon catcher air gap. NYI
void SetKalFwdAir(float kalfwdair)
distance (dist,cm) through air traveled, when projected forwards to detector edge. For ND only, muon catcher air gap. NYI
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
void SetTrkBakCell(int trkbakcell)
void SetVtxDist(float vtxdist)
minimum distance of start position of best Kalman track to detector wall. Negative indicates outside ...
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
void SetKalBakCell(int kalbakcell)
note that fwd<->bak values are interchanged if track direction is backwards
void SetCosFwdCellND(int cosfwdcellnd)
cells in that line segment instead of distance of line segment, including MC, ND only
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
float getActivity(art::PtrVector< rb::CellHit > const &trackHits, art::PtrVector< rb::CellHit > const &sliceHits)
void SetKalFwdCellND(int kalfwdcellnd)
cells in that line segment instead of distance of line segment, including MC, ND only
double ProjectedAirDist(TVector3 vertex, TVector3 dir)
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
void SetFScattMax(float fscattmax)
maximum scattering angle (Fernanda's)
void SetEndDist(float enddist)
minimum distance of end position of best Kalman track to detector wall. Negative indicates outside of...
void SetCosHitRatio(float coshitratio)
ratio of hits in cosmic track to hits in slice
art::ServiceHandle< cheat::BackTracker > bt
void SetAngleCos(float anglel)
cos of angle of Cosmic track wrt beam
void SetTrkBakSteel(float trkbaksteel)
distance (dist,cm) through air traveled, when projected forwards to detector edge. For ND only, muon catcher air gap. NYI
const art::ProductToken< std::vector< rb::Cluster > > fSlicerToken
virtual void InterpolateXY(double z, double &x, double &y) const
void SetTrkBakDist(float trkbakdist)
distance (dist,cm) of best track from end to det wall projected forwards (Fwd), based on Kalman Track...
TMVA::Reader * fReaderTAP1
cosrej::CosRejFxs * fCosRejFxs
TMVA::Reader * fReaderUCTuned
TMVA::Reader * fReaderTAP2
T min(const caf::Proxy< T > &a, T b)
bool IsNoise() const
Is the noise flag set?
TVector3 Stop() const
Position of the final trajectory point.
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
void SetCosChisq(float coschisq)
same for Cosmic track
void SetTrkEPerNHit(float trkepernhit)
Track energy per kalman track hit.
void SetCosSlope(float cosslope)
same for Cosmic track
virtual double DistanceFromEnd(double z) const
void SetKalBakDist(float kalbakdist)
cosmic track distance projected backwards from track start
ProductToken< T > consumes(InputTag const &)
void SetCosBakDist(float cosbakdist)
cosmic track distance projected backwards from track start
void SetFScattSig(float fscattsig)
sigma of scattering angles (Fernanda's);
void SetAngleKal2nd(float angle2nd)
cos of angle of 2nd best Kalman track wrt beam
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
std::string fGeneratorLabel
unsigned int HighestPIDTrack(const std::vector< art::Ptr< rb::Track > > &sliceTracks, const std::string &remidModuleLabel, const art::Event &e)
art::ServiceHandle< geo::Geometry > geom
void SetConCosPID(float concospid)
PID for contained cosmic rejection (The one to use - P5 train)
void SetFScattExt(float fscattext)
scattering variable (Fernanda's)
bool fTrackContainmentOnly
void SetTrkFwdDist(float trkfwddist)
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.
void SetERatio(float eratio)
ratio of visible E in best track/ visible E in whole slice
void SetKalFwdSteel(float kalfwdsteel)
distance (dist,cm) through muon catcher traveled, when projected forwards to detector edge...
void SetCosBakSteel(float cosbaksteel)
distance (dist,cm) through muon catcher traveled, when projected backwards to detector edge...
TMVA::Reader * fReaderProd5_Per1