19 #include "NovaDAQConventions/DAQConventions.h" 21 #include "Utilities/func/MathUtil.h" 37 fPIDModuleLabel (pset.
get<
std::
string>(
"PIDModuleLabel" )),
38 fPIDMethod (pset.
get< bool >(
"PIDMethod" )),
39 fTMVAName (pset.
get<
std::
string>(
"TMVAFileName",
"UC_Eres_BDTG.weights.xml")),
40 fTMVANameSingle (pset.
get<
std::
string>(
"TMVAFileSingle",
"UC_MuonE_Single_BDTG.weights.xml")),
41 fTMVANameNonSingle (pset.
get<
std::
string>(
"TMVAFileNonSingle",
"UC_MuonE_NonSingle_BDTG.weights.xml")),
86 fReaderE->BookMVA(
"BDTG",fileName.c_str());
250 double trackLen = 0.0;
255 if(bestTrack == -1)
return sliceEnergy;
256 if(bestTrack == 999)
return sliceEnergy;
258 if (
fPIDMethod) trackLen = sliceTracks[bestTrack]->TotalLength();
261 double recoTrkGeV = sliceTracks[bestTrack]->TotalGeV();
274 trackHits = sliceTracks[bestTrack]->AllCells();
280 double vertexEnergy = 0.0;
281 if (vertexHits.
NCell()>0){
282 vertexEnergy = vertexHits.
TotalGeV();
286 double extraHadE = trkCleanUpAlg->
ExtraEOnTrackInGeV(*sliceTracks[bestTrack],*sliceCluster);
289 double totHadE = vertexEnergy + extraHadE;
290 double hadQEGeV =
HadEQE(totHadE);
292 double hadCCGeV =
HadECC(totHadE,
fRun, ismc, isRHC);
296 double calCCMuonWeight = 1.81134;
297 double calCCVertWeight = 2.17923;
298 double calcc = calCCMuonWeight*recoTrkGeV + calCCVertWeight*vertexEnergy;
304 double qeE = muonGeVFromLen + hadQEGeV;
307 double nonqeE = muonGeVFromLen + hadNonQEGeV;
310 double ccE = muonGeVFromLen + hadCCGeV;
315 sliceEnergy.
SetE(ccE);
346 double trackLen = 0.0;
353 if(bestTrack == -1)
return 0;
354 if(bestTrack == 999)
return 0;
356 if (
fPIDMethod) trackLen = sliceTracks[bestTrack]->TotalLength();
360 TVector3 mudir = sliceTracks[bestTrack]->Dir();
361 if(mudir.Mag()!=1) mudir = mudir.Unit();
380 double Mn = 0.939565;
381 double Mnp = Mn - Mb;
382 double pmu_2 = muE*muE - (Mu*Mu);
384 double recoNu = (2*Mnp*muE-(Mnp*Mnp+Mu*Mu-(Mp*Mp))) / (2*(Mnp-muE+
sqrt(pmu_2)*costheta));
387 double errMb = 0.005;
390 if (trackType ==
kFDMuon) {resE = 0.043;}
397 double errE = resE*muE;
399 double rescostheta = 0.05;
402 double varMb = errMb*errMb;
403 double varE = errE*errE;
404 double varcostheta = (rescostheta*costheta)*(rescostheta*costheta);
406 double x = 2*Mnp*muE - (Mnp*Mnp + Mu*Mu - (Mp*Mp));
407 double varx = 4*(Mn*Mn+Mb*Mb)*varE + 4*(Mn*Mn+muE*muE +0.5*Mb*Mb)*varMb;
408 double vary = 4*varMb + 4*varE + 4*((muE*muE*costheta*costheta*varE)/(2*pmu_2) + pmu_2*varcostheta);
410 double sigmaE_2 = recoNu*recoNu*((varx+recoNu*recoNu*vary)/(x*x));
411 error =
sqrt(sigmaE_2);
488 double muonWeight_act = 0.900903;
489 double muonWeight_tran = 3.84225;
490 double muonWeight_steel = 6.58836;
491 double vertWeight = 1.61335;
492 double offsetWeight = 0.985003;
497 double trackLen = 0.0;
505 if(bestTrack == -1)
return sliceEnergy;
506 if(bestTrack == 999)
return sliceEnergy;
508 if (
fPIDMethod) trackLen = sliceTracks[bestTrack]->TotalLength();
513 double recoTrkGeV_act = 0.0;
514 double recoTrkGeV_tran = 0.0;
515 double recoTrkGeV_steel = 0.0;
516 double actTrkLen = 0.0;
517 double catTrkLen = 0.0;
518 double tranX = 999.0;
519 double tranY = 999.0;
520 double hadCalGeV_act = 0.0;
521 double hadCalGeV_tran = 0.0;
522 double hadCalGeV_steel = 0.0;
523 double hadTrkGeV_act = 0.0;
524 double hadTrkGeV_tran = 0.0;
525 double hadTrkGeV_steel = 0.0;
526 double vertexEnergy = 0.0;
527 double vertexTrkEnergy = 0.0;
533 NDTrackLength(sliceTracks[bestTrack],&actTrkLen,&catTrkLen,&tranX,&tranY);
536 trackHits = sliceTracks[bestTrack]->AllCells();
542 if (vertexHits.
NCell()>0){
543 vertexEnergy = vertexHits.
TotalGeV();
547 for (
unsigned int hitHad = 0; hitHad <vertexHits.
NCell(); ++hitHad){
550 hadCalGeV_act += vertexHits.
RecoHit(hitHad).
GeV();
553 hadCalGeV_tran += vertexHits.
RecoHit(hitHad).
GeV();
556 hadCalGeV_steel += vertexHits.
RecoHit(hitHad).
GeV();
564 for (std::map<int,float>::iterator
it = extraHadE.begin();
it!=extraHadE.end(); ++
it){
566 hadTrkGeV_act +=
it->second;
569 hadTrkGeV_tran +=
it->second;
572 hadTrkGeV_steel +=
it->second;
576 vertexTrkEnergy = hadTrkGeV_act + hadTrkGeV_tran + hadTrkGeV_steel;
588 double totHadE = vertexEnergy + vertexTrkEnergy;
589 double hadQEGeV =
NDHadEQE(totHadE);
591 double hadCCGeV =
NDHadECC(totHadE, isRHC);
595 double calcc = muonWeight_act*recoTrkGeV_act + muonWeight_tran*recoTrkGeV_tran + muonWeight_steel*recoTrkGeV_steel
596 + vertWeight*vertexEnergy + offsetWeight;
599 double qeE = muonGeVFromLen + hadQEGeV;
602 double nonqeE = muonGeVFromLen + hadNonQEGeV;
605 double ccE = muonGeVFromLen + hadCCGeV;
610 sliceEnergy.
SetE(ccE);
646 if(
fND)
return NDEnergy(sliceTracks, sliceHits, sliceCluster, e, trkCleanUpAlg, isRHC);
647 if(
fFD)
return FDEnergy(sliceTracks, sliceHits, sliceCluster, e, trkCleanUpAlg, isRHC);
657 const std::vector<rb::CellHit>& allHits,
669 double trackLen = -1;
677 if(bestTrack == -1)
return truthEnergy;
678 if(bestTrack == 999)
return truthEnergy;
686 if (funcReturn.empty())
return truthEnergy;
689 std::set<int> neutrinoTrackIDs;
690 for (
unsigned int i = 0;
i < particleList.size(); ++
i){
691 neutrinoTrackIDs.insert(particleList[
i]->TrackId());
694 std::set<int>::iterator itr = neutrinoTrackIDs.begin();
696 bool goodMuon =
false;
697 double trueCatcherE = -5.0;
698 double trueMuonE = -5.0;
700 while( itr != neutrinoTrackIDs.end() ){
706 if ((PDG != 13)&&(PDG != -13)){
713 if (mother->
TrackId() != id_int){
726 trueMuonE = particle->
E();
730 double closeToCenter = 50.0;
734 for(
unsigned int q = 0; q < flshits.size(); ++q){
736 if(
fabs(flshits[q].GetEntryY()) < closeToCenter){
738 trueCatcherE = flshits[q].GetEntryEnergy() + 0.105658;
739 closeToCenter =
fabs(flshits[q].GetEntryY());
741 if(
fabs(flshits[q].GetExitY()) < closeToCenter){
742 trueCatcherE = flshits[q].GetExitEnergy() + 0.105658;
743 closeToCenter =
fabs(flshits[q].GetExitY());
768 for (
unsigned int sliceHit = 0; sliceHit < sliceHits.
size(); ++sliceHit){
772 for (
unsigned int trackHit = 0; trackHit < trackHits.
size(); ++trackHit){
773 match = (sliceHits[sliceHit] == trackHits[trackHit]);
777 if (!match) vertexCluster.
Add(sliceHits[sliceHit]);
781 return vertexCluster;
787 double& longestTrackLength)
const 791 for(
size_t c = 0;
c < sliceTracks.size(); ++
c){
800 if ((sliceTracks[
c]->TotalLength() > longestTrackLength)&&(sliceTracks[
c]->TotalLength() <
fMaxTrkLen)) {
801 longestTrackLength = sliceTracks[
c]->TotalLength();
816 double* transitionY)
const 819 std::cout<<
"You have a problem! Your track only has one trajectory point!"<<
std::endl;
823 bool foundActRegion =
false;
824 bool foundCatRegion =
false;
826 double activeLen = 0.0;
827 double catcherLen = 0.0;
832 double tranX = 999.0;
833 double tranY = 999.0;
837 for(
int n = 0;
n < N-1; ++
n){
840 foundActRegion =
true;
844 foundCatRegion =
true;
846 foundActRegion =
true;
860 if (foundActRegion && foundCatRegion){
861 TVector3
diff = catTraj-actTraj;
862 double addActLen = (
fTranZPos-actTraj.Z())*diff.Mag()/diff.Z();
863 double addCatLen = (catTraj.Z()-
fTranZPos)*diff.Mag()/diff.Z();
864 activeLen += addActLen;
865 catcherLen += addCatLen;
868 tranX = actTraj.X() + (
fTranZPos - actTraj.Z())*diff.X()/diff.Z();
869 tranY = actTraj.Y() + (
fTranZPos - actTraj.Z())*diff.Y()/diff.Z();
872 if (actLength) *actLength = activeLen;
873 if (catLength) *catLength = catcherLen;
874 if (transitionX) *transitionX = tranX;
875 if (transitionY) *transitionY = tranY;
883 double* trkCalCat)
const 885 double trkEnergyInActive = 0.0;
886 double trkEnergyInTran = 0.0;
887 double trkEnergyInCatcher = 0.0;
890 for (
unsigned int hitIdx = 0; hitIdx <track->
NCell(); ++hitIdx){
893 trkEnergyInActive += track->
RecoHit(hitIdx).
GeV();
896 trkEnergyInTran += track->
RecoHit(hitIdx).
GeV();
899 trkEnergyInCatcher += track->
RecoHit(hitIdx).
GeV();
904 if (trkCalAct) *trkCalAct = trkEnergyInActive;
905 if (trkCalTran) *trkCalTran = trkEnergyInTran;
906 if (trkCalCat) *trkCalCat = trkEnergyInCatcher;
919 if (trkLen <= 0.0)
return muonE;
963 if (muonE <= 0.0)
return trkLen;
999 if ((catTrkLen <= 0.0) && (actTrkLen <= 0.0))
return muonE;
1018 if (((trkCalAct > 0.0)||(trkCalTran > 0.0))&&(trkCalCat == 0.0)){
1023 if ((trkCalAct == 0.0)&&(trkCalTran == 0.0)&&(trkCalCat > 0.0)){
1028 if (((trkCalAct > 0.0)||(trkCalTran > 0.0))&&(trkCalCat > 0.0)){
1048 double actTrkLen = 0.0;
1049 double catTrkLen = 0.0;
1050 double trkCalAct = 0.0;
1051 double trkCalTran = 0.0;
1052 double trkCalCat = 0.0;
1087 if (catTrkLen <= 0.0)
return muonE;
1089 double slope = 0.00513341;
1090 double offset = 0.230294;
1092 muonE = slope*catTrkLen +
offset;
1110 if (catTrkLen <= 0.0)
return muonE;
1148 std::vector<rb::Energy> energiesOfTracks;
1150 for(
size_t c = 0;
c < sliceTracks.size(); ++
c){
1158 energiesOfTracks.push_back(trackEnergy);
1165 double trkLen = sliceTracks[
c]->TotalLength();
1167 trackEnergy.
SetE(energy);
1168 energiesOfTracks.push_back(trackEnergy);
1173 trackEnergy.
SetE(energy);
1174 energiesOfTracks.push_back(trackEnergy);
1179 return energiesOfTracks;
1212 if (calE < 0.0)
return hadE;
1214 double slope1 = 0.621;
1215 double slope2 = 1.47;
1216 double slope3 = 1.81;
1217 double slope4 = 2.06;
1219 double stitch1 = 0.0597;
1220 double stitch2 = 0.139;
1221 double stitch3 = 0.800;
1223 if (calE < stitch1){
1224 hadE = slope1*calE +
offset;
1226 else if (calE < stitch2){
1227 hadE = slope2*calE + ((slope1-slope2)*stitch1 + offset);
1229 else if (calE < stitch3){
1230 hadE = slope3*calE + ((slope2-slope3)*stitch2 + (slope1-slope2)*stitch1 +
offset);
1233 hadE = slope4*calE + ((slope3-slope4)*stitch3 + (slope2-slope3)*stitch2
1234 +(slope1-slope2)*stitch1 + offset);
1250 if (calE < 0.0)
return hadE;
1252 double slope1 = 0.94;
1253 double slope2 = 1.21;
1254 double slope3 = 1.78;
1255 double slope4 = 2.20;
1257 double stitch1 = 0.132;
1258 double stitch2 = 0.223;
1259 double stitch3 = 0.964;
1261 if (calE < stitch1){
1262 hadE = slope1*calE +
offset;
1264 else if (calE < stitch2){
1265 hadE = slope2*calE + ((slope1-slope2)*stitch1 + offset);
1267 else if (calE < stitch3){
1268 hadE = slope3*calE + ((slope2-slope3)*stitch2 + (slope1-slope2)*stitch1 +
offset);
1271 hadE = slope4*calE + ((slope3-slope4)*stitch3 + (slope2-slope3)*stitch2
1272 +(slope1-slope2)*stitch1 + offset);
1296 if (calE < 0.0)
return hadE;
1298 double slope1 = 1.567;
1299 double slope2 = 1.822;
1301 double stitch1 = 0.0820;
1303 if (calE < stitch1){
1304 hadE = slope1*calE +
offset;
1307 hadE = slope2*calE + ((slope1-slope2)*stitch1 + offset);
1324 if (calE < 0.0)
return hadE;
1326 double slope1 = 1.080;
1327 double slope2 = 1.889;
1329 double stitch1 = 0.169;
1331 if (calE < stitch1){
1332 hadE = slope1*calE +
offset;
1335 hadE = slope2*calE + ((slope1-slope2)*stitch1 + offset);
1353 std::vector< art::Ptr<cosrej::CosRejObj> > cosmicrej = CRO.at(slicenum);
1354 Evars[0] = totalgev;
1355 Evars[1] = cosmicrej[0]->AngleKal();
1356 Evars[2] = cosmicrej[0]->Scatt();
1357 Evars[3] = cosmicrej[0]->ERatio();
1358 Evars[4] = cosmicrej[0]->KalFwdCell();
1359 Evars[5] = cosmicrej[0]->CosHitRatio()*1.0*nhitslice;
1360 Evars[6] = cosmicrej[0]->NTracks3D();
1361 return fReaderE->EvaluateRegression(
"BDTG")[0];
1372 double trackLen = 0.0;
1373 double hadclustcalE = 0.0;
1378 if(bestTrack == -1)
return -5.0;
1379 if(bestTrack == 999)
return -5.0;
1382 trackHits = sliceTracks[bestTrack]->AllCells();
1387 if (!(vertexHits.
NCell()==0)){
1394 std::vector< art::Ptr<cosrej::CosRejObj> > cosmicrej = CRO.at(slicenum);
1398 TMVAvarsSingle[3] = (cosmicrej[0]->Scatt())/(cosmicrej[0]->TrackLenKal());;
1411 double trackLen = 0.0;
1412 double hadclustcalE = 0.0;
1417 if(bestTrack == -1)
return -5.0;
1418 if(bestTrack == 999)
return -5.0;
1421 trackHits = sliceTracks[bestTrack]->AllCells();
1426 if (!(vertexHits.
NCell()==0)){
1433 std::vector< art::Ptr<cosrej::CosRejObj> > cosmicrej = CRO.at(slicenum);
1448 std::vector< art::Ptr<cosrej::CosRejObj> > cosmicrej = CRO.at(slicenum);
1450 pars[0] = cosmicrej[0]->FScattExt();
1451 pars[1] = cosmicrej[0]->FScattSig();
1452 pars[2] = cosmicrej[0]->FScattMax();
1453 pars[3] = cosmicrej[0]->FScattSum();
1455 double scatterMuEnergy = nnEnergy->
Value(0, pars);
1457 return scatterMuEnergy;
double E(const int i=0) const
size_t NTrajectoryPoints() const
double fTranZPos
Z position of middle of transition plane.
void NDTrackEnergySplitting(art::Ptr< rb::Track > const &track, double *trkCalAct=NULL, double *trkCalTran=NULL, double *trkCalCat=NULL) const
Function for ND that divides track visible energy into three regions: active, transiton, and catcher.
double GetUCMuonENonSingle(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, art::Event const &e, art::FindManyP< cosrej::CosRejObj > CRO, size_t slicenum)
Get values from cosrej object and use them to get TMVA trained value of uncontained muon energy for F...
back track the reconstruction to the simulation
double NDMuonCatcherEForActiveAndCatcher(double catTrkLen) const
Function that, given muon track length in cm, returns muon energy in GeV. This is only for ND muons t...
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
void GetCenter(double *xyz, double localz=0.0) const
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Energy estimators for CC events.
TMVA::Reader * fReaderUCMuonSingle
Reader for TMVA to get muon uncontained energy.
double FDMslope3
GeV/cm; Third slope for 4 spline fit for muon trk len -> energy in FD.
fvar< T > fabs(const fvar< T > &x)
double GetFSE(art::FindManyP< cosrej::CosRejObj > CRO, size_t slicenum)
Fernanda's scattering energy. Scattering variables themselves stored in CosRej.
double FDMstitch3
cm; Third stitch location on 4 spline fit for muon trk len -> energy in FD
float Evars[8]
Inputs for uncontained energy BDT (TMVA)
double FDMslope2
GeV/cm; Second slope for 4 spline fit for muon trk len -> energy in FD.
const sim::Particle * TrackIDToMotherParticle(int const &id) const
double HadENonQE(double calE) const
Function that, given total hadron cal energy on and off muon track in GeV, returns weighted hadron en...
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...
double NDHadECC(double calE, bool isRHC) const
Function that, given total hadron cal energy on and off muon track in GeV, returns weighted hadron en...
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
double predict_prod5_fd_muon_energy(double trkLen, int run, bool ismc, bool isRHC)
unsigned short Plane() const
void SetNDTrkCalAct(float ndtrkcalact)
double NDHadEQE(double calE) const
Function that, given total hadron cal energy on and off muon track in GeV, returns weighted hadron en...
void SetMCGoodTrueMuon(bool mcgoodtruemuon)
std::string fTMVANameNonSingle
path to TMVA muon weight file
NumuE Energy(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, art::Ptr< rb::Cluster > sliceCluster, art::Event const &e, murem::TrackCleanUpAlg *trkCleanUpAlg, bool isRHC)
Returns various energy estimations for the current detector. See comments at top of function in ...
double FDMuonEFromTrackLength(double trkLen, int run, bool ismc, bool isRHC) const
Function that, given muon track length in cm, returns muon energy in GeV. This is done with a spline ...
double predict_prod5_nd_had_energy(double hadVisE, bool isRHC)
const CellGeo * Cell(int icell) const
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
NumuE FDEnergy(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, art::Ptr< rb::Cluster > sliceCluster, art::Event const &e, murem::TrackCleanUpAlg *trkCleanUpAlg, bool isRHC)
Returns various energy estimations for FD. See comments at top of function in .cxx for full explanati...
void SetMCTrueMuonE(float mctruemuone)
void SetTrkNonQEE(float trknonqe)
A collection of associated CellHits.
int LongestTrack(std::vector< art::Ptr< rb::Track > > const &sliceTracks, double &longestTrackLength) const
Function that finds the longest track that is still smaller than maxTrkLen.
NumuEAlg(fhicl::ParameterSet const &pset)
int NDTransitionPlane() const
Function that returns the plane number of the transition plane - the last active plane before the ste...
double predict_prod5_nd_cat_energy(double ndtrklencat, bool isRHC)
void SetRecoTrkCCHadE(float recotrkcchade)
const PlaneGeo * Plane(unsigned int i) const
double HadECC(double calE, int run, bool ismc, bool isRHC) const
Function that, given total hadron cal energy on and off muon track in GeV, returns weighted hadron en...
double NDMCoffset
GeV; Offset for linear fit for muon trk len -> energy in ND muon catcher.
std::string find_file(std::string const &filename) const
double NDMuonEInActiveAndCatcher(double actTrkLen, double catTrkLen, bool isRHC) const
Function that, given muon track length in cm, returns muon energy in GeV. This is only for ND muons t...
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
double MuonEFromTrackLength(double trkLen) const
Function that, given muon track length in cm, returns muon energy in GeV. This is done with a spline ...
int isnan(const stan::math::var &a)
void NDTrackLength(art::Ptr< rb::Track > const &track, double *actLength=NULL, double *catLength=NULL, double *transitionX=NULL, double *transitionY=NULL) const
Function for ND that divides track length into length in active material and length in muon catcher...
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
double FDMoffset
GeV; Offset for 4 spline fit for muon trk len -> energy in FD.
void SetNDHadCalCat(float ndhadcalcat)
void SetMCTrueMuonCatcherE(float mctruemuoncatchere)
std::map< int, float > ExtraEOnTrackPlanesInGeV(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
double NDMCslope
GeV/cm; First slope for linear fit for muon trk len -> energy in ND muon catcher. ...
Far Detector at Ash River, MN.
void SetNDHadCalAct(float ndhadcalact)
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
int fTranPlane
Last plane of scintillator before steel begins in ND.
void SetNDHadTrkCat(float ndhadtrkcat)
TVector3 fNumiDir
Direction of NuMI Beam.
double GetUCMuonESingle(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, art::Event const &e, art::FindManyP< cosrej::CosRejObj > CRO, size_t slicenum)
Get values from cosrej object and use them to get TMVA trained value of uncontained muon energy for F...
double NDMuonEFromTrackLength(double actTrkLen, double catTrkLen, double trkCalAct, double trkCalTran, double trkCalCat, bool isRHC, MuonType *muonType=NULL) const
Function that, given muon track length in cm, returns muon energy in GeV. This distinguishes between ...
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
void SetTrkQEE(float trkqe)
void SetNDHadTrkAct(float ndhadtrkact)
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Near Detector in the NuMI cavern.
double FDMslope1
GeV/cm; First slope for 4 spline fit for muon trk len -> energy in FD.
rb::Energy QEFormulaEnergy(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, double &error, art::Event const &e, bool isRHC)
Returns QE energy estimation using formula. See comments at top of function in .cxx for full explanat...
void SetNDTrkTranX(float ndtrktranx)
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
std::string fPIDModuleLabel
Label for module writing ReMId objects to file.
TMVA::Reader * fReaderE
Reader for TMVA to get uncontained energy.
NumuE MCTruthEnergyVariables(const std::vector< art::Ptr< rb::Track >> &sliceTracks, const art::Ptr< rb::Cluster > &slice, const std::vector< rb::CellHit > &allHits, const art::Event &e)
Returns various useful truth energy values, good for fitting.
void SetTrkCCE(float trkcce)
double DetHalfHeight() const
double NDTransitionZPos() const
Function that returns the z position of the middle of the transition plane - the last active plane be...
bool checkIsMC(art::Event const &evt) const
void SetNDTrkLenCat(float ndtrklencat)
void SetNDTrkCalTran(float ndtrkcaltran)
bool fPIDMethod
If true, use highest PID to identify muon track instead of longest track (Only FD currently) ...
double FDMstitch2
cm; Second stitch location on 4 spline fit for muon trk len -> energy in FD
double ActiveTrackLengthFromMuonE(double muonE) const
Function that, given muon energy in GeV, returns an effective track length in active material in cm...
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
std::string fTMVANameSingle
path to TMVA muon weight file
double GetUCE(art::FindManyP< cosrej::CosRejObj > CRO, size_t slicenum, float totalgev, int nhitslice)
Get values from cosrej object and use them to get TMVA trained value of uncontained energy for FD...
TMVA::Reader * fReaderUCMuonNonSingle
Reader for TMVA to get muon uncontained energy.
void SetNDHadTrkTran(float ndhadtrktran)
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
void SetNDTrkTranY(float ndtrktrany)
rb::Cluster MakeVertexCluster(art::PtrVector< rb::CellHit > const &sliceHits, art::PtrVector< rb::CellHit > const &trackHits) const
Function that makes cluster from all hits in slice that do not belong to track.
void SetNDTrkLenAct(float ndtrklenact)
double DetHalfWidth() const
A container for energy information.
double pythag(double x, double y)
2D Euclidean distance
double FDMslope4
GeV/cm; Fourth slope for 4 spline fit for muon trk len -> energy in FD.
void SetHadTrkE(float hadtrk)
double NDMuonEInActiveOnly(double actTrkLen, bool isRHC) const
Function that, given muon track length in cm, returns muon energy in GeV. This is only for ND muons t...
void SetCalCCE(float calcc)
void SetNDHadCalTran(float ndhadcaltran)
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
NumuE NDEnergy(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, art::Ptr< rb::Cluster > sliceCluster, art::Event const &e, murem::TrackCleanUpAlg *trkCleanUpAlg, bool isRHC)
Returns various energy estimations for ND. See comments at top of function in .cxx for full explanati...
void SetHadCluster(rb::Cluster hadcluster)
double fMaxTrkLen
Maximum track length allowed in detector (diagonal length)
bool PassMuonCuts(int trackID, art::PtrVector< rb::CellHit > const &hits) const
Tool for NumuEAna that allows one to see if primary muon (or any track ID you feed the function) cont...
double HadEQE(double calE) const
Function that, given total hadron cal energy on and off muon track in GeV, returns weighted hadron en...
double predict_prod5_fd_had_energy(double hadVisE, int run, bool ismc, bool isRHC)
void SetRecoMuonE(float recomuone)
float TMVAvarsSingle[38]
Inputs for muon uncontained energy BDT (TMVA) for single track events.
void SetHadCalE(float hadcal)
void SetNDTrkCalCat(float ndtrkcalcat)
float TMVAvarsNonSingle[38]
Inputs for muon uncontained energy BDT (TMVA) for non-single track events.
double predict_prod5_nd_act_energy(double ndtrklenact, bool isRHC)
double FDMstitch1
cm; First stitch location on 4 spline fit for muon trk len -> energy in FD
float ExtraEOnTrackInGeV(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
std::string fTMVAName
path to TMVA weight file
double Value(int index, double in0, double in1, double in2, double in3)
double NDHadENonQE(double calE) const
Function that, given total hadron cal energy on and off muon track in GeV, returns weighted hadron en...
Encapsulate the geometry of one entire detector (near, far, ndos)
unsigned int HighestPIDTrack(const std::vector< art::Ptr< rb::Track > > &sliceTracks, const std::string &remidModuleLabel, const art::Event &e)
std::vector< rb::Energy > MuonEnergies(std::vector< art::Ptr< rb::Track > > const sliceTracks, int run, bool ismc, bool isRHC) const
Function that returns energy for each track given the assumption that it is a muon. It uses the functions above, like MuonEFromTrackLength and NDMuonEFromTrackLength.
double NDMuonEInCatcherOnly(double catTrkLen) const
Function that, given muon track length in cm, returns muon energy in GeV. This is only for ND muons t...
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.