Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | List of all members
fnex::NuMuAnalysisSetup Class Reference

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-12-01/FNEX/custom/NuMuAnalysisSetup.h"

Inheritance diagram for fnex::NuMuAnalysisSetup:
fnex::AnalysisSetupBase

Public Member Functions

 NuMuAnalysisSetup (fhicl::ParameterSet const &pset)
 
virtual ~NuMuAnalysisSetup ()=default
 
void FillVars (art::Event const &evt, art::ValidHandle< std::vector< rb::Cluster > > const &slices, art::InputTag const &tag, std::vector< fnex::DataVarVals > &dataVarVals, std::vector< fnex::MCVarVals > &mcVarVals, std::vector< fnex::EventId > &evIds, fnex::BeamType_t const &beamType) const override
 

Protected Member Functions

bool FillTruthVars (fnex::MCVarVals &vars, art::Event const &evt, size_t const &sliceNum, art::FindOneP< simb::MCTruth > const &fomct, art::FindOne< simb::MCFlux > const &fomcf, art::FindOne< fxwgt::FluxWeights > const &fofw, art::FindOne< simb::GTruth > const &fogt) const
 
void FillGENIEWeightVars (fnex::MCVarVals &vars, rwgt::GENIEReweightTable const &) const
 
double FindXSecWeight (simb::MCTruth const &nu, simb::GTruth const &gt, rwgt::WeightType const &tag) const
 
double FindXSecWeightSA (simb::MCTruth const &nu) const
 
void FindRPAWeights (fnex::MCVarVals &vars, simb::MCTruth const &mct, simb::GTruth const &gt, rwgt::WeightType const &tag) const
 
void FindMECWeight (fnex::MCVarVals &vars, simb::MCTruth const &mct, simb::GTruth const &gt, rwgt::WeightType const &tag) const
 
void FindDISWeight (fnex::MCVarVals &vars, simb::MCTruth const &mct, simb::GTruth const &gt, rwgt::WeightType const &tag) const
 

Protected Attributes

art::ServiceHandle< art::TFileServicefTFS
 TFileService. More...
 
art::ServiceHandle< rwgt::MCReweightfRwgt
 MCReweight service. More...
 
rwgt::WeightType fXSecWgtScheme
 which "Tufts" weighting to use More...
 

Private Member Functions

void Reconfigure (fhicl::ParameterSet const &pset) override
 
void FillRecoVars (fnex::DataVarVals &vars, numue::NumuE const &, rb::Energy const &, rb::Track const &, bool isMC, fnex::BeamType_t const &beamType) const
 
void SecondAnalysisEnergy (numue::NumuE const &numue, rb::Track const &track, double &hadE, double &muonE) const
 
void ThirdAnalysisEnergy (numue::NumuE const &numue, rb::Track const &track, double &hadE, double &muonE, bool isMC, fnex::BeamType_t const &beamType) const
 
void Analysis2018Energy (numue::NumuE const &numue, rb::Track const &track, double &hadE, double &muonE, bool isMC, fnex::BeamType_t const &beamType) const
 
int ThirdAnalysisQuantile (double &recoE, double &hadEFrac) const
 
int FourthAnalysisQuantile (double &lepE, double &hadE, fnex::BeamType_t const &beamType) const
 
int FourthAnalysisQuantileFHC (double &lepE, double &hadE) const
 
int FourthAnalysisQuantileRHC (double &lepE, double &hadE) const
 
unsigned int HighestPIDTrack (std::vector< art::Ptr< rb::Track > > const &sliceTracks, art::InputTag const &remidTag, art::Event const &e) const
 

Private Attributes

std::string fNuMuEEstimator
 
int fHadEFracQuantile
 Hadronic Energy Fraction quantile. More...
 
art::ServiceHandle< nova::dbi::RunHistoryServicefRH
 RunHistory service. More...
 
art::ServiceHandle< ds::DetectorServicefDetService
 which detector? More...
 

Detailed Description

Definition at line 36 of file NuMuAnalysisSetup.h.

Constructor & Destructor Documentation

fnex::NuMuAnalysisSetup::NuMuAnalysisSetup ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 36 of file NuMuAnalysisSetup.cxx.

References Reconfigure().

37  {
38  this->Reconfigure(pset);
39  }
void Reconfigure(fhicl::ParameterSet const &pset) override
virtual fnex::NuMuAnalysisSetup::~NuMuAnalysisSetup ( )
virtualdefault

Member Function Documentation

void fnex::NuMuAnalysisSetup::Analysis2018Energy ( numue::NumuE const &  numue,
rb::Track const &  track,
double &  hadE,
double &  muonE,
bool  isMC,
fnex::BeamType_t const &  beamType 
) const
private

Definition at line 922 of file NuMuAnalysisSetup.cxx.

References fillBadChanDBTables::det, ds::DetectorService::DetId(), fDetService, fRH, numue::NumuE::HadCalE(), numue::NumuE::HadTrkE(), make_mec_shifts_plots::isRHC, novadaq::cnv::kFARDET, novadaq::cnv::kNEARDET, fnex::kRHC, LOG_DEBUG, numue::NumuE::NDTrkLenAct(), numue::NumuE::NDTrkLenCat(), NumuEnergyFunc::predict_prod4_fd_had_energy(), NumuEnergyFunc::predict_prod4_fd_muon_energy(), NumuEnergyFunc::predict_prod4_nd_act_energy(), NumuEnergyFunc::predict_prod4_nd_cat_energy(), NumuEnergyFunc::predict_prod4_nd_had_energy(), nova::dbi::RunHistory::RunNumber(), and rb::Track::TotalLength().

Referenced by FillRecoVars().

928  {
929  auto det = fDetService->DetId();
930  auto run = fRH->RunNumber();
931 
932  double hadVisE = numue.HadTrkE() + numue.HadCalE();
933 
934  bool isRHC = false; // FHC mode
935  if(beamType == fnex::kRHC) isRHC = true ; // RHC mode
936  LOG_DEBUG("NuMuAnalysisSetup")
937  << "Inside NuMuAnalysisSetup::Analysis2018Energy function. ";
938  //if(isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "RHC file";
939  //else LOG_VERBATIM("NuMuAnalysisSetup") << "FHC file";
940 
942  {
943 
944  double trkLen = track.TotalLength();
945 
946  //if(isMC && !isRHC) LOG_DEBUG("NuMuAnalysisSetup") << "MC, not RHC";
947  //if(!isMC && !isRHC) LOG_DEBUG("NuMuAnalysisSetup") << "Not MC, not RHC";
948 
949  muonE = NumuEnergyFunc::predict_prod4_fd_muon_energy(trkLen, run, isMC, isRHC);
950  hadE = NumuEnergyFunc::predict_prod4_fd_had_energy(hadVisE, run, isMC, isRHC);
951 
952  LOG_DEBUG("NuMuAnalysisSetup")
953  << "Inside NuMuAnalysisSetup::Analysis2018Energy function. "
954  << "LepE = " << muonE << "\t HadE = " << hadE;
955  }
956 
958  {
959 
960  double trkLenAct = numue.NDTrkLenAct();
961  double trkLenCat = numue.NDTrkLenCat();
962 
963  // if(isMC && !isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "MC, not RHC";
964  // if(!isMC && !isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "Not MC, not RHC";
965 
967  hadE = NumuEnergyFunc::predict_prod4_nd_had_energy(hadVisE, isRHC);
968  }
969 
970  }// end Ana2018 energy estimator
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
novadaq::cnv::DetId DetId() const
What detector are we in?
Energy estimators for CC events.
Definition: FillEnergies.h:7
Definition: event.h:19
Far Detector at Ash River, MN.
art::ServiceHandle< nova::dbi::RunHistoryService > fRH
RunHistory service.
Near Detector in the NuMI cavern.
double predict_prod4_nd_cat_energy(double ndtrklencat, bool isRHC)
double predict_prod4_nd_act_energy(double ndtrklenact, bool isRHC)
Definition: run.py:1
double predict_prod4_fd_muon_energy(double trkLen, int run, bool ismc, bool isRHC)
Definition: Numu2018Fits.cxx:7
art::ServiceHandle< ds::DetectorService > fDetService
which detector?
double predict_prod4_fd_had_energy(double hadVisE, int run, bool ismc, bool isRHC)
double predict_prod4_nd_had_energy(double hadVisE, bool isRHC)
int RunNumber() const
Definition: RunHistory.h:374
void fnex::AnalysisSetupBase::FillGENIEWeightVars ( fnex::MCVarVals vars,
rwgt::GENIEReweightTable const &  rwgtTable 
) const
protectedinherited

Definition at line 47 of file AnalysisSetupBase.cxx.

References fnex::AnalysisSetupBase::fXSecWgtScheme, fnex::MCVarVals::GENIEWgts(), fnex::isSumSmallXSecJoint2017Name(), fnex::isUnusedGENIEName(), rwgt::kXSecCVWgt_2017, LOG_DEBUG, LOG_WARNING, rwgt::GENIEReweightTable::Minus1Sigma(), rwgt::GENIEReweightTable::Minus2Sigma(), rwgt::GENIEReweightTable::NShifts(), rwgt::GENIEReweightTable::Plus1Sigma(), rwgt::GENIEReweightTable::Plus2Sigma(), fnex::MCVarVals::set_GENIEWgt(), and string.

Referenced by fnex::AnalysisSetupBase::FillTruthVars().

49  {
50  static std::map<rwgt::EReweightLabel, std::string> labelMap
51  {
52  // NC elastic parameters
53  {rwgt::fReweightMaNCEL, "MaNCEL"},
54  {rwgt::fReweightEtaNCEL, "EtaNCEL"},
55 
56  // CCQE parameters
57  {rwgt::fReweightNormCCQE, "NormCCQE"},
58  {rwgt::fReweightNormCCQEenu, "NormCCQE_EnuDep"},
59  {rwgt::fReweightMaCCQEshape, "MaCCQE_shape"},
60  {rwgt::fReweightMaCCQE, "MaCCQE"},
61  {rwgt::fReweightVecCCQEshape, "VecCCQE_shape"},
62 
63  // Resonant production: CC
64  {rwgt::fReweightNormCCRES, "NormCCRES"},
65  {rwgt::fReweightMaCCRESshape, "MaCCRES_shape"},
66  {rwgt::fReweightMvCCRESshape, "MvCCRES_shape"},
67  {rwgt::fReweightMaCCRES, "MaCCRES"},
68  {rwgt::fReweightMvCCRES, "MvCCRES"},
69 
70  // Resonant production: NC
71  {rwgt::fReweightNormNCRES, "NormNCRES"},
72  {rwgt::fReweightMaNCRESshape, "MaNCRES_shape"},
73  {rwgt::fReweightMvNCRESshape, "MvNCRES_shape"},
74  {rwgt::fReweightMaNCRES, "MaNCRES"},
75  {rwgt::fReweightMvNCRES, "MvNCRES"},
76 
77  // Coherent pion production
78  {rwgt::fReweightMaCOHpi, "MaCOHpi"},
79  {rwgt::fReweightR0COHpi, "R0COHpi"},
80 
81  // Non-resonant background
82  {rwgt::fReweightRvpCC1pi, "RvpCC1pi"},
83  {rwgt::fReweightRvpCC2pi, "RvpCC2pi"},
84  {rwgt::fReweightRvpNC1pi, "RvpNC1pi"},
85  {rwgt::fReweightRvpNC2pi, "RvpNC2pi"},
86  {rwgt::fReweightRvnCC1pi, "RvnCC1pi"},
87  {rwgt::fReweightRvnCC2pi, "RvnCC2pi"},
88  {rwgt::fReweightRvnNC1pi, "RvnNC1pi"},
89  {rwgt::fReweightRvnNC2pi, "RvnNC2pi"},
90  {rwgt::fReweightRvbarpCC1pi, "RvbarpCC1pi"},
91  {rwgt::fReweightRvbarpCC2pi, "RvbarpCC2pi"},
92  {rwgt::fReweightRvbarpNC1pi, "RvbarpNC1pi"},
93  {rwgt::fReweightRvbarpNC2pi, "RvbarpNC2pi"},
94  {rwgt::fReweightRvbarnCC1pi, "RvbarnCC1pi"},
95  {rwgt::fReweightRvbarnCC2pi, "RvbarnCC2pi"},
96  {rwgt::fReweightRvbarnNC1pi, "RvbarnNC1pi"},
97  {rwgt::fReweightRvbarnNC2pi, "RvbarnNC2pi"},
98 
99  // DIS tweaking parameters
100  {rwgt::fReweightAhtBY, "AhtBY"},
101  {rwgt::fReweightBhtBY, "BhtBY"},
102  {rwgt::fReweightCV1uBY, "CV1uBY"},
103  {rwgt::fReweightCV2uBY, "CV2uBY"},
104  {rwgt::fReweightAhtBYshape, "AhtBYshape"},
105  {rwgt::fReweightBhtBYshape, "BhtBYshape"},
106  {rwgt::fReweightCV1uBYshape, "CV1uBYshape"},
107  {rwgt::fReweightCV2uBYshape, "CV2uBYshape"},
108  {rwgt::fReweightNormDISCC, "NormDISCC"},
109  {rwgt::fReweightRnubarnuCC, "RnubarnuCC"},
110  {rwgt::fReweightDISNuclMod, "DISNuclMod"},
111 
112  {rwgt::fReweightNC, "NC"},
113 
114  // Hadronization (free nucleon target)
115  {rwgt::fReweightAGKY_xF1pi, "AGKY_xF1pi"},
116  {rwgt::fReweightAGKY_pT1pi, "AGKY_pT1pi"},
117 
118  // Medium-effects to hadronization
119  {rwgt::fReweightFormZone, "FormZone"},
120 
121  // Intranuclear rescattering parameters.
122  {rwgt::fReweightMFP_pi, "MFP_pi"},
123  {rwgt::fReweightMFP_N, "MFP_N"},
124  {rwgt::fReweightFrCEx_pi, "FrCEx_pi"},
125  {rwgt::fReweightFrInel_pi, "FrInel_pi"},
126  {rwgt::fReweightFrAbs_pi, "FrAbs_pi"},
127  {rwgt::fReweightFrPiProd_pi, "FrPiProd_pi"},
128  {rwgt::fReweightFrCEx_N, "FrCEx_N"},
129  // the next two appear no longer to be in use with GENIE v3
130  //{rwgt::fReweightFrElas_pi, "FrElas_pi"},
131  //{rwgt::fReweightFrElas_N, "FrElas_N"},
132  {rwgt::fReweightFrInel_N, "FrInel_N"},
133  {rwgt::fReweightFrAbs_N, "FrAbs_N"},
134  {rwgt::fReweightFrPiProd_N, "FrPiProd_N"},
135 
136  // Nuclear model
137  {rwgt::fReweightCCQEPauliSupViaKF, "CCQEPauliSupViaKF"},
138  {rwgt::fReweightCCQEMomDistroFGtoSF, "CCQEMomDistroFGtoSF"},
139 
140  // Resonance decays
141  {rwgt::fReweightBR1gamma, "BR1gamma"},
142  {rwgt::fReweightBR1eta, "BR1eta"},
143  {rwgt::fReweightTheta_Delta2Npi, "Theta_Delta2Npi"}
144  };
145 
146  std::string rwgtLabel;
147  for (const auto rwgtPair : labelMap){
148  unsigned int rwgtInt = rwgtPair.first; // note: casting from enum to int since the rwgtTable wants ints
149  rwgtLabel = rwgtPair.second;
150  //To do: switch this to log_debug when you're confident this all works
151  LOG_DEBUG("AnalysisSetupBase")
152  << "Looking at GENIE var " << rwgtLabel
153  << " (" << rwgtTable.Minus2Sigma(rwgtInt) << ", "
154  << rwgtTable.Minus1Sigma(rwgtInt) << ", "
155  << rwgtTable.Plus1Sigma(rwgtInt) << ", "
156  << rwgtTable.Plus2Sigma(rwgtInt) << ")";
157  if (static_cast<unsigned int>(rwgtInt) > rwgtTable.NShifts())
158  continue;
159 
161  rwgtLabel = "SumSmallXSecJoint2017";
162 
163  // don't bother with GENIE weights we don't currently use
164  if( fnex::isUnusedGENIEName(rwgtLabel) ){
165  LOG_DEBUG("AnalysisSetupBase") << "\tactually don't store it, it's not used ";
166  continue;
167  }
168 
169  LOG_DEBUG("AnalysisSetupBase") << "\tstore as " << rwgtLabel;
170 
171  // check that the weight for this variable and event isn't completely crazy.
172  // It is, don't add the weight to the vars
173  if(rwgtTable.Minus2Sigma(rwgtInt) > 100. ||
174  rwgtTable.Minus1Sigma(rwgtInt) > 100. ||
175  rwgtTable.Plus1Sigma (rwgtInt) > 100. ||
176  rwgtTable.Plus2Sigma (rwgtInt) > 100. ){
177 
178  LOG_WARNING("AnalysisSetupBase")
179  << "Weight for "
180  << rwgtLabel
181  << " is crazy for at least one shift, skip this knob for this event";
182 
183  continue;
184  }
185 
186  // the set_GENIEWgt function multiplies the existing values already in
187  // the GENIEVar by the table vals for each sigma so when we are setting
188  // the SumSmallXSecJoint2017 variable it becomes an product of weights from all the
189  // individual contributions. Also, only the SumSmallXSecJoint2017 key will be stored
190  // in the eventual EventList tree file.
191 
192  vars.set_GENIEWgt(rwgtLabel, "-2sigma", rwgtTable.Minus2Sigma(rwgtInt));
193  vars.set_GENIEWgt(rwgtLabel, "-1sigma", rwgtTable.Minus1Sigma(rwgtInt));
194  vars.set_GENIEWgt(rwgtLabel, "+1sigma", rwgtTable.Plus1Sigma (rwgtInt));
195  vars.set_GENIEWgt(rwgtLabel, "+2sigma", rwgtTable.Plus2Sigma (rwgtInt));
196 
197  LOG_DEBUG("AnalysisSetupBase")
198  << "GENIE var " << rwgtLabel
199  << " stored as ("
200  << vars.GENIEWgts(rwgtLabel).at(-2.) << ", "
201  << vars.GENIEWgts(rwgtLabel).at(-1.) << ", "
202  << vars.GENIEWgts(rwgtLabel).at( 1.) << ", "
203  << vars.GENIEWgts(rwgtLabel).at( 2.) << ")";
204 
205  }
206 
207  return;
208  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
bool isUnusedGENIEName(std::string const &varname)
Definition: VarVals.h:238
std::map< float, float > GENIEWgts(unsigned char const &varkey)
Definition: VarVals.cxx:116
rwgt::WeightType fXSecWgtScheme
which "Tufts" weighting to use
#define LOG_WARNING(category)
void set_GENIEWgt(unsigned char const &varkey, std::string const &sigma, float const &wgt)
Definition: VarVals.cxx:282
bool isSumSmallXSecJoint2017Name(std::string const &varname)
Definition: VarVals.h:193
enum BeamMode string
void fnex::NuMuAnalysisSetup::FillRecoVars ( fnex::DataVarVals vars,
numue::NumuE const &  numue,
rb::Energy const &  energy,
rb::Track const &  track,
bool  isMC,
fnex::BeamType_t const &  beamType 
) const
private

Definition at line 170 of file NuMuAnalysisSetup.cxx.

References Analysis2018Energy(), fillBadChanDBTables::det, ds::DetectorService::DetId(), rb::Prong::Dir(), rb::Energy::E(), e, fDetService, fNuMuEEstimator, make_mec_shifts_plots::isRHC, novadaq::cnv::kFARDET, fnex::kHad_RecoE, fnex::kLep_RecoE, fnex::kLep_RecoE_MCFrac, novadaq::cnv::kNEARDET, fnex::kNu_RecoE, fnex::kNuMuSelection, fnex::kRecoQ2, fnex::kRHC, LOG_DEBUG, fetch_tb_beamline_files::md, numue::NumuE::NDTrkLenAct(), numue::NumuE::NDTrkLenCat(), cet::pow(), predict_prod3_nd_act_energy(), predict_prod3_nd_cat_energy(), SecondAnalysisEnergy(), fnex::MetaData::selectionType, fnex::DataVarVals::set_val_at(), std::sqrt(), ThirdAnalysisEnergy(), numue::NumuE::TrkCCE(), and fnex::DataVarVals::val_at().

Referenced by FillVars().

176  {
177  // lepton energy is muon energy
178  double hadE = (numue.TrkCCE() - energy.E());
179  double muonE = energy.E();
180  if(!fNuMuEEstimator.compare("SA"))
181  this->SecondAnalysisEnergy(numue, track, hadE, muonE);
182  else if(!fNuMuEEstimator.compare("TA"))
183  this->ThirdAnalysisEnergy(numue, track, hadE, muonE, isMC, beamType);
184  else if(!fNuMuEEstimator.compare("2018")){
185  LOG_DEBUG("NuMuAnalysisSetup") << "Using 2018 numu energy estimator";
186  this->Analysis2018Energy(numue, track, hadE, muonE, isMC, beamType);
187  }
188  else
189  throw cet::exception("NuMuAnalysisSetup:FillRecoVars : Need to set NuMuEEstimator to SA or TA (for 2nd and 3rd analysis respectively)");
190 
191  vars.set_val_at(fnex::kLep_RecoE, muonE);
192  vars.set_val_at(fnex::kHad_RecoE, hadE );
193 
194  // fraction of leptonic energy in muon catcher for ND
195  auto det = fDetService->DetId();
196  //auto run = fRH->RunNumber();
197 
198  if(det == novadaq::cnv::kNEARDET) {
199 
200  bool isRHC = false; // FHC mode
201  if(beamType == fnex::kRHC) isRHC = true ; // RHC mode
202 
203  double lepEAct = predict_prod3_nd_act_energy(numue.NDTrkLenAct(), isRHC);
204  double lepECat = predict_prod3_nd_cat_energy(numue.NDTrkLenCat(), isRHC);
205  double lepMCFrac = lepECat / (lepEAct + lepECat);
206  vars.set_val_at(fnex::kLep_RecoE_MCFrac, lepMCFrac);
207  }
208  else vars.set_val_at(fnex::kLep_RecoE_MCFrac, 0);
209 
210  // reconstructed Q^2
211  double M_mu_sqr = std::pow(0.1056, 2);
212  double p_mu = std::sqrt(std::pow(energy.E(), 2) - M_mu_sqr);
213 
214  TVector3 beamDir;
216  beamDir = TVector3(-8.424e-04,
217  -6.174e-02,
218  +9.981e-1).Unit();
219  else if(det == novadaq::cnv::kFARDET)
220  beamDir = TVector3(-6.833e-05,
221  +6.388e-02,
222  +9.980e-1).Unit();
223 
226 
227  double RecoQ2 = 2 * vars.val_at(fnex::kNu_RecoE, md) * (energy.E() - p_mu * track.Dir().Dot(beamDir) - M_mu_sqr);
228 
229  vars.set_val_at(fnex::kRecoQ2, RecoQ2);
230 
231  return;
232  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
novadaq::cnv::DetId DetId() const
What detector are we in?
static const unsigned char kNu_RecoE
Definition: VarVals.h:36
Energy estimators for CC events.
Definition: FillEnergies.h:7
void ThirdAnalysisEnergy(numue::NumuE const &numue, rb::Track const &track, double &hadE, double &muonE, bool isMC, fnex::BeamType_t const &beamType) const
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
static const unsigned char kRecoQ2
Definition: VarVals.h:38
Definition: event.h:19
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void Analysis2018Energy(numue::NumuE const &numue, rb::Track const &track, double &hadE, double &muonE, bool isMC, fnex::BeamType_t const &beamType) const
Far Detector at Ash River, MN.
static const unsigned char kLep_RecoE_MCFrac
Definition: VarVals.h:37
double energy
Definition: plottest35.C:25
Near Detector in the NuMI cavern.
float val_at(unsigned char const &varkey, fnex::MetaData const &md) const
Definition: VarVals.cxx:499
art::ServiceHandle< ds::DetectorService > fDetService
which detector?
fnex::SelectionType_t selectionType
Definition: Structs.h:52
double predict_prod3_nd_cat_energy(double ndtrklencat, bool isRHC)
void set_val_at(std::string const &varkey, float const &val)
Definition: VarVals.cxx:525
static const unsigned char kHad_RecoE
Definition: VarVals.h:34
Float_t e
Definition: plot.C:35
static const unsigned char kLep_RecoE
Definition: VarVals.h:35
void SecondAnalysisEnergy(numue::NumuE const &numue, rb::Track const &track, double &hadE, double &muonE) const
double predict_prod3_nd_act_energy(double ndtrklenact, bool isRHC)
bool fnex::AnalysisSetupBase::FillTruthVars ( fnex::MCVarVals vars,
art::Event const &  evt,
size_t const &  sliceNum,
art::FindOneP< simb::MCTruth > const &  fomct,
art::FindOne< simb::MCFlux > const &  fomcf,
art::FindOne< fxwgt::FluxWeights > const &  fofw,
art::FindOne< simb::GTruth > const &  fogt 
) const
protectedinherited

Definition at line 211 of file AnalysisSetupBase.cxx.

References evt, fnex::AnalysisSetupBase::FillGENIEWeightVars(), fnex::AnalysisSetupBase::FindDISWeight(), fnex::AnalysisSetupBase::FindMECWeight(), fnex::AnalysisSetupBase::FindRPAWeights(), fnex::AnalysisSetupBase::FindXSecWeight(), simb::MCFlux::fndecay, simb::MCFlux::fntype, simb::GTruth::fNumPi0, simb::GTruth::fNumPiMinus, simb::GTruth::fNumPiPlus, simb::MCFlux::fptype, simb::MCFlux::ftptype, simb::MCFlux::ftpx, simb::MCFlux::ftpy, simb::MCFlux::ftpz, fnex::AnalysisSetupBase::fXSecWgtScheme, simb::MCTruth::GetNeutrino(), fxwgt::FluxWeights::GetTotalCentralValue(), gt(), art::Event::isRealData(), cet::maybe_ref< T >::isValid(), simb::kCC, simb::kDIS, simb::kElectronScattering, simb::kMEC, simb::kRes, fnex::kTrueCCNC, fnex::kTrueE, fnex::kTrueHitNuc, fnex::kTrueIntMode, fnex::kTrueIntType, fnex::kTrueParentDecay, fnex::kTrueParentPDG, fnex::kTrueParentPT, fnex::kTrueParentPZ, fnex::kTrueParentTargetPDG, fnex::kTruePDG, fnex::kTruePDGOrig, fnex::kXSecCVPPFX_Weight, LOG_DEBUG, LOG_WARNING, gen_flatrecord::pt, cet::maybe_ref< T >::ref(), fnex::MCVarVals::set_val_at(), std::sqrt(), and fnex::MCVarVals::val_at().

Referenced by FillVars(), and fnex::NuEAnalysisSetup::FillVars().

218  {
219  // do nothing if this is real data, return true because we did what we
220  // needed to
221  if( evt.isRealData() ) return true;
222 
223  if(fomct.size() < sliceNum){
224  LOG_WARNING("AnalysisSetupBase")
225  << "The association between slices and MCTruth objects "
226  << "is too small for the requested slice number "
227  << fomct.size()
228  << " / "
229  << sliceNum
230  << " bail on this slice";
231  return false;
232  }
233 
234  art::Ptr<simb::MCTruth> truth = fomct.at(sliceNum);
235 
236  LOG_DEBUG("AnalysisSetupBase")
237  << "Current MCTruth is "
238  << *truth;
239 
240  // if we can't find an MCTruth object we are in trouble, so move on
241  if(!truth){
242  LOG_WARNING("AnalysisSetupBase")
243  << "No MCTruth found for slice number "
244  << sliceNum
245  << " bail on this slice";
246  return false;
247  }
248 
249  // now get the FluxWeights and GTruth
250  double fluxWeight = 1.;
251  cet::maybe_ref<fxwgt::FluxWeights const> fw = fofw.at(sliceNum);
252  if(!fw.isValid()){
253  LOG_WARNING("AnalysisSetupBase")
254  << "No FluxWeights found for slice number "
255  << sliceNum;
256  }
257  else{
258  fluxWeight = fw.ref().GetTotalCentralValue();
259  LOG_DEBUG("AnalysisSetupBase")
260  << "FluxWeight set to "
261  << fluxWeight;
262  }
263 
264  cet::maybe_ref<simb::GTruth const> gt = fogt.at(sliceNum);
265  if(!gt.isValid())
266  LOG_WARNING("AnalysisSetupBase")
267  << "No GTruth found for slice number "
268  << sliceNum;
269 
270 
271  // Fill the information related to the MCTruth object
272  auto nu = truth->GetNeutrino();
273 
274  vars.set_val_at(fnex::kTrueE, nu.Nu().E() );
275  vars.set_val_at(fnex::kTruePDG, nu.Nu().PdgCode() );
276  vars.set_val_at(fnex::kTrueCCNC, nu.CCNC() );
277  vars.set_val_at(fnex::kTrueIntType, nu.InteractionType());
278  vars.set_val_at(fnex::kTrueIntMode, nu.Mode() );
279 
280  // HitNuc is used for the MECInitStateNPFracWeight systematic
281  vars.set_val_at(fnex::kTrueHitNuc, nu.HitNuc() );
282 
283  // get the GENIEReweight informaiton
285 
286  // For a numu or nutau scattering off an electron, the exchanged particle
287  // is a Z. But for a nue, there are two diagrams, one with a Z and one
288  // with a W. Both the final states have the electron ejected, so the
289  // interaction is a quantum-mechanical superposition between the
290  // two. It's neither really CC or NC. We do need to set this boolean
291  // though. Treating nu-on-e scattering as NC (ie not oscillating it) is
292  // wrong, the cross-section for nue is different to numu and nutau so
293  // oscillations do have an effect on the final spectrum. If we mark it CC
294  // the different oscillated fluxes will be weighted together in the fit
295  // appropriately, so this is safe, even if not perfectly physically
296  // accurate.
297  //
298  // Above reason is given by CAF group in the CAFMaker_module.cc
299  // We are simply taking this argument to make things consistent between FNEX and CAFAna.
301 
302  LOG_DEBUG("AnalysisSetupBase")
303  << " Inside AnalysisSetupBase::FillTruthVars function "
304  << " Run : "
305  << evt.run()
306  << " Subrun : "
307  << evt.subRun()
308  << " Event : "
309  << evt.id().event()
310  << " Slice : "
311  << sliceNum
312  << " TrueE "
313  << nu.Nu().E()
314  << " tru->GetNeutrino().CCNC() : "
315  << nu.CCNC()
316  << " srTruth.iscc : "
317  << vars.val_at(fnex::kTrueCCNC)
318  << " tru->GetNeutrino().Mode() : "
319  << nu.Mode()
320  << " kElectronScattering : "
322  << " simb::kCC : "
323  << simb::kCC;
324 
325  // In the 2017 version, the XsecCV weight is the product of the RPA,
326  // EmpericalMEC and non-res 1Pi. We'll also put the flux weights into
327  // the kXSecCVPPFX_Weight value for the event
328 
329  if( frw.isValid() && frw.size() > 0){
330  art::Ptr<rwgt::GENIEReweightTable> rw = frw.at(0); // should only be one...
331  if(rw) this->FillGENIEWeightVars(vars, *rw);
332 
334  this->FindXSecWeight(*truth, gt.ref(), fXSecWgtScheme) * fluxWeight);
335 
336  LOG_DEBUG("AnalysisSetupBase")
337  << "Check Weights "
338  << this->FindXSecWeight(*truth, gt.ref(), fXSecWgtScheme)
339  << " "
340  << fluxWeight
341  << " "
343 
344  LOG_DEBUG("AnalysisSetupBase")
345  << "About to start loading up 2017 systematics!"
346  << "\n \t CC? " << vars.val_at(fnex::kTrueCCNC) << " == " << 1. * simb::kCC
347  << "\n \t RES? " << vars.val_at(fnex::kTrueIntType) << " == " << 1. * simb::kRes
348  << "\n \t DIS? " << vars.val_at(fnex::kTrueIntType) << " == " << 1. * simb::kDIS
349  << "\n \t MEC? " << vars.val_at(fnex::kTrueIntType) << " == " << 1. * simb::kMEC
350  << "\n \t nu/nubar? " << vars.val_at(fnex::kTruePDG) << " > 0"
351  << "\n \t npions = " << gt.ref().fNumPiPlus + gt.ref().fNumPiMinus + gt.ref().fNumPi0
352  << "\n \t hits a neutron? " << nu.HitNuc() << " == 2112";
353 
354  this->FindRPAWeights(vars, *truth, gt.ref(), fXSecWgtScheme);
355  this->FindMECWeight (vars, *truth, gt.ref(), fXSecWgtScheme);
356  this->FindDISWeight (vars, *truth, gt.ref(), fXSecWgtScheme);
357 
358  }
359  else
360  LOG_WARNING("AnalysisSetupBase")
361  << "No GENIEReweightTable found for slice number "
362  << sliceNum;
363 
364  // now get the flux information
365  if(fomcf.size() < sliceNum)
366  LOG_WARNING("AnalysisSetupBase")
367  << "The association between MCFlux and slices is "
368  << "too small for the requested slice number "
369  << fomcf.size()
370  << " / "
371  << sliceNum
372  << " skip Flux information for this slice";
373  else{
374  cet::maybe_ref<simb::MCFlux const> mcf = fomcf.at(sliceNum);
375  if( mcf.isValid() ){
376  float pt = std::sqrt(mcf.ref().ftpx * mcf.ref().ftpx +
377  mcf.ref().ftpy * mcf.ref().ftpy);
378 
379  vars.set_val_at(fnex::kTruePDGOrig, mcf.ref().fntype );
383  vars.set_val_at(fnex::kTrueParentPZ, mcf.ref().ftpz );
384  vars.set_val_at(fnex::kTrueParentPT, pt );
385  }
386  else
387  LOG_WARNING("AnalysisSetupBase")
388  << "No MCFlux found for slice number "
389  << sliceNum;
390  } // end if fomcf is too small
391 
392  return true;
393  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
double FindXSecWeight(simb::MCTruth const &nu, simb::GTruth const &gt, rwgt::WeightType const &tag) const
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
double ftpx
Definition: MCFlux.h:79
static const unsigned char kTrueParentPT
Definition: VarVals.h:48
int ftptype
Definition: MCFlux.h:82
void set_val_at(std::string const &varkey, float const &val)
Definition: VarVals.cxx:164
T sqrt(T number)
Definition: d0nt_math.hpp:156
static const unsigned char kXSecCVPPFX_Weight
Definition: VarVals.h:56
void FindRPAWeights(fnex::MCVarVals &vars, simb::MCTruth const &mct, simb::GTruth const &gt, rwgt::WeightType const &tag) const
static const unsigned char kTrueHitNuc
Definition: VarVals.h:52
int fNumPiPlus
number of pi pluses after reaction, before FSI
Definition: GTruth.h:78
int fNumPiMinus
number of pi minuses after reaction, before FSI
Definition: GTruth.h:79
static const unsigned char kTrueParentPDG
Definition: VarVals.h:47
bool isValid() const
Definition: maybe_ref.h:60
static const unsigned char kTruePDGOrig
Definition: VarVals.h:46
void FindMECWeight(fnex::MCVarVals &vars, simb::MCTruth const &mct, simb::GTruth const &gt, rwgt::WeightType const &tag) const
static const unsigned char kTrueCCNC
Definition: VarVals.h:44
rwgt::WeightType fXSecWgtScheme
which "Tufts" weighting to use
int fptype
Definition: MCFlux.h:63
int fndecay
Definition: MCFlux.h:50
int fNumPi0
number of pi0 after reaction, before FSI
Definition: GTruth.h:77
int evt
float val_at(unsigned char const &varkey) const
Access a Var by name.
Definition: VarVals.cxx:77
static const unsigned char kTrueIntType
Definition: VarVals.h:45
double ftpz
Definition: MCFlux.h:81
static const unsigned char kTrueParentDecay
Definition: VarVals.h:50
#define LOG_WARNING(category)
static const unsigned char kTrueParentTargetPDG
Definition: VarVals.h:51
int fntype
Definition: MCFlux.h:51
static const unsigned char kTrueE
Definition: VarVals.h:42
void FillGENIEWeightVars(fnex::MCVarVals &vars, rwgt::GENIEReweightTable const &) const
static const unsigned char kTrueIntMode
Definition: VarVals.h:53
double ftpy
Definition: MCFlux.h:80
bool gt(unsigned short int a, unsigned short int b)
static const unsigned char kTrueParentPZ
Definition: VarVals.h:49
void FindDISWeight(fnex::MCVarVals &vars, simb::MCTruth const &mct, simb::GTruth const &gt, rwgt::WeightType const &tag) const
float GetTotalCentralValue() const
Get the total central value correction calculated as the product of the HP and Beam Focusing componen...
Definition: FluxWeights.cxx:46
static const unsigned char kTruePDG
Definition: VarVals.h:43
void fnex::NuMuAnalysisSetup::FillVars ( art::Event const &  evt,
art::ValidHandle< std::vector< rb::Cluster > > const &  slices,
art::InputTag const &  tag,
std::vector< fnex::DataVarVals > &  dataVarVals,
std::vector< fnex::MCVarVals > &  mcVarVals,
std::vector< fnex::EventId > &  evIds,
fnex::BeamType_t const &  beamType 
) const
overridevirtual

Implements fnex::AnalysisSetupBase.

Definition at line 42 of file NuMuAnalysisSetup.cxx.

References art::Event::event(), fHadEFracQuantile, FillRecoVars(), fnex::AnalysisSetupBase::FillTruthVars(), PandAna.Demos.pi0_spectra::fmt, FourthAnalysisQuantile(), art::FindManyP< ProdB, Data >::get(), remid::HighestPIDTrack(), art::Event::isRealData(), cet::maybe_ref< T >::isValid(), fnex::kHad_RecoE, fnex::kLep_RecoE, LOG_DEBUG, LOG_WARNING, cet::sqlite::max(), min(), cet::maybe_ref< T >::ref(), art::Event::run(), datagram_client::sl, art::Event::subRun(), and fnex::DataVarVals::val_at().

49  {
50  dataVarVals.clear();
51  mcVarVals .clear();
52  evIds .clear();
53 
54  if( !slices.isValid() ){
55  LOG_WARNING("NuMuAnalysisSetup")
56  << "could not find valid handle of rb::Cluster, ignore this event";
57  }
58  else{
59 
60  LOG_DEBUG("NuMuAnalysisSetup")
61  << "Start filling variables for tag: "
62  << tag
63  << " and quantile "
65 
66  // get the FindOnes
67  art::FindOne <numue::NumuE > fone (slices, evt, tag);
69  art::FindOneP <simb::MCTruth > fomct(slices, evt, tag);
70  art::FindOne <simb::MCFlux > fomcf(slices, evt, tag);
72  art::FindOne <simb::GTruth > fogt (slices, evt, tag);
73 
74  std::vector< art::Ptr<rb::Track> > trackPtrs;
75 
76  size_t bestPIDTrack = std::numeric_limits<size_t>::max();
77 
78  for(size_t sl = 0; sl < slices->size(); ++sl){
79  fnex::DataVarVals dataVars;
80  fnex::MCVarVals mcVars;
81  fnex::EventId evId(evt.run(),
82  evt.subRun(),
83  evt.event(),
85  (*slices)[sl].ID());
86 
87  // check that the FindOnes are valid
88  // only need the 0th entry in the FindOnes because we have only one slice here
89  if(fone.isValid() ){
90 
91  cet::maybe_ref<numue::NumuE const> nume = fone.at(sl);
92 
93  fmt.get(sl, trackPtrs);
94  if(trackPtrs.size() < 1) continue;
95 
96  // get the best pid for the slice - this index will not be the
97  // same as in the preskimmed slice
98  bestPIDTrack = remid::HighestPIDTrack(trackPtrs, tag, evt);
99  //bestPIDTrack = this->HighestPIDTrack(trackPtrs, tag, evt); // Reimplement here because Is3D does not work
100  // and the remid default does not exclude 2d tracks for Skimmed files
101 
102  // this is not terribly efficient, as we are making the FindOneP
103  // for every slice in the event, but I don't see a good way around
104  // that while still getting the bookkeeping right
105  art::FindOne<rb::Energy> foe(trackPtrs, evt, tag);
106 
107  if( nume.isValid() && foe.isValid() )
108  this->FillRecoVars(dataVars, nume.ref(), foe.at(bestPIDTrack).ref(), *(trackPtrs[bestPIDTrack]), !evt.isRealData(), beamType);
109 
110  }// end filling recovariables
111  else{
112  LOG_WARNING("NuMuAnalysisSetup")
113  << "something is wrong with the FindOne<numue::NumuE> "
114  << " skip this slice";
115  continue;
116  }
117 
118  // now for the MC variables, if this is a MC event
119  if( !this->FillTruthVars(mcVars, evt, sl, fomct, fomcf, fofw, fogt) ) continue;
120 
121  LOG_DEBUG("NuMuAnalysisSetup")
122  << "numu truth is "
123  << mcVars;
124 
125  // check the hadronic eenrgy fraction quantile and only add the vars to the collection
126  // if we are in the desired range (expects 0 for SA
127  double lepE = dataVars.val_at(fnex::kLep_RecoE, fnex::MetaData());
128  double hadE = dataVars.val_at(fnex::kHad_RecoE, fnex::MetaData());
129 
130  int hadEFracQ = FourthAnalysisQuantile(lepE, hadE, beamType);
131  if( fHadEFracQuantile == 0 || hadEFracQ == fHadEFracQuantile ){
132  dataVarVals.push_back(dataVars);
133  mcVarVals .push_back(mcVars);
134  evIds .push_back(evId);
135  }
136  } // end loop over slices
137 
138  } // end if the handle was valid
139 
140  return;
141  } // NuMuAnalysisSetup::FillVars
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
int fHadEFracQuantile
Hadronic Energy Fraction quantile.
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:930
bool isValid() const
Definition: maybe_ref.h:60
int evt
void FillRecoVars(fnex::DataVarVals &vars, numue::NumuE const &, rb::Energy const &, rb::Track const &, bool isMC, fnex::BeamType_t const &beamType) const
float val_at(unsigned char const &varkey, fnex::MetaData const &md) const
Definition: VarVals.cxx:499
int FourthAnalysisQuantile(double &lepE, double &hadE, fnex::BeamType_t const &beamType) const
#define LOG_WARNING(category)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
static const unsigned char kHad_RecoE
Definition: VarVals.h:34
bool FillTruthVars(fnex::MCVarVals &vars, art::Event const &evt, size_t const &sliceNum, art::FindOneP< simb::MCTruth > const &fomct, art::FindOne< simb::MCFlux > const &fomcf, art::FindOne< fxwgt::FluxWeights > const &fofw, art::FindOne< simb::GTruth > const &fogt) const
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
static const unsigned char kLep_RecoE
Definition: VarVals.h:35
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:874
unsigned int HighestPIDTrack(const std::vector< art::Ptr< rb::Track > > &sliceTracks, const std::string &remidModuleLabel, const art::Event &e)
Definition: ReMId.cxx:249
bool isValid() const
Definition: Handle.h:326
void fnex::AnalysisSetupBase::FindDISWeight ( fnex::MCVarVals vars,
simb::MCTruth const &  mct,
simb::GTruth const &  gt,
rwgt::WeightType const &  tag 
) const
protectedinherited

Definition at line 604 of file AnalysisSetupBase.cxx.

References simb::GTruth::fNumPi0, simb::GTruth::fNumPiMinus, simb::GTruth::fNumPiPlus, simb::MCTruth::GetNeutrino(), simb::kCCDIS, fnex::kDISvnCC1pi_Weight, simb::kNCDIS, fnex::kOtherDIS_Weight, fnex::kTrueIntType, fnex::kTruePDG, rwgt::kXSecCVWgt_2017, rwgt::kXSecCVWgt_2018, fnex::MCVarVals::set_val_at(), and fnex::MCVarVals::val_at().

Referenced by fnex::AnalysisSetupBase::FillTruthVars().

608  {
609  auto nu = mct.GetNeutrino();
610 
611  // 2017 DISvnCC1pi
612  // 2017 DIS weights: DISvnCC1pi gets its own weight
613  // All other v/vbar n/p NC/CC 0-4pi combinations get lumped into SumSmallXSecNumu2017
614  vars.set_val_at(fnex::kDISvnCC1pi_Weight, 1.); // will be changed if event is DIS CC neutrino hitting neutron producing 1 pion
615  vars.set_val_at(fnex::kOtherDIS_Weight, 1.); // will be changed if event is DIS but doesn't fall into the above category
616 
617  if(tag == rwgt::kXSecCVWgt_2017 ||
619 
620  auto numPions = gt.fNumPiPlus + gt.fNumPiMinus + gt.fNumPi0;
621 
622  if(vars.val_at(fnex::kTrueIntType) == 1. * simb::kCCDIS && //CC and DIS
623  vars.val_at(fnex::kTruePDG) > 0 && //not anti-nu
624  numPions == 1 && //1 pion
625  nu.HitNuc() == 2112){ //hits a neutron
626  if(nu.W() < 3.) vars.set_val_at(fnex::kDISvnCC1pi_Weight, 1.5);
627  else vars.set_val_at(fnex::kDISvnCC1pi_Weight, 1.05);
628  } else if((vars.val_at(fnex::kTrueIntType) == 1. * simb::kCCDIS ||
629  vars.val_at(fnex::kTrueIntType) == 1. * simb::kNCDIS ) && //DIS
630  (nu.HitNuc() == 2112 ||
631  nu.HitNuc() == 2212 ) && //hits a neutron or a proton
632  numPions < 4){ // <4 pions
633  if(nu.W() < 3.) vars.set_val_at(fnex::kOtherDIS_Weight, 1.5);
634  else vars.set_val_at(fnex::kOtherDIS_Weight, 1.05);
635  }
636  }
637 
638  return;
639  }
charged current deep inelastic scatter
Definition: MCNeutrino.h:134
static const unsigned char kOtherDIS_Weight
Definition: VarVals.h:63
static const unsigned char kDISvnCC1pi_Weight
Definition: VarVals.h:59
void set_val_at(std::string const &varkey, float const &val)
Definition: VarVals.cxx:164
float val_at(unsigned char const &varkey) const
Access a Var by name.
Definition: VarVals.cxx:77
static const unsigned char kTrueIntType
Definition: VarVals.h:45
charged current deep inelastic scatter
Definition: MCNeutrino.h:133
bool gt(unsigned short int a, unsigned short int b)
static const unsigned char kTruePDG
Definition: VarVals.h:43
void fnex::AnalysisSetupBase::FindMECWeight ( fnex::MCVarVals vars,
simb::MCTruth const &  mct,
simb::GTruth const &  gt,
rwgt::WeightType const &  tag 
) const
protectedinherited

Definition at line 545 of file AnalysisSetupBase.cxx.

References fnex::AnalysisSetupBase::fRwgt, simb::MCTruth::GetNeutrino(), rwgt::MCReweight::GetWeight(), rwgt::kEmpiricalMEC_2017, rwgt::kEmpiricalMEC_2018, fnex::kEmpiricalMEC_Weight, fnex::kEmpiricalMECtoGENIEQE_Weight, fnex::kEmpiricalMECtoGENIERES_Weight, simb::kMEC, fnex::kTruePDG, rwgt::kXSecCVWgt_2017, rwgt::kXSecCVWgt_2018, LOG_DEBUG, rwgt::MCReweight::MECq0q3ShapeSyst2018(), fnex::MCVarVals::set_val_at(), and fnex::MCVarVals::val_at().

Referenced by fnex::AnalysisSetupBase::FillTruthVars().

549  {
550  auto nu = mct.GetNeutrino();
551 
552  // 2017 MEC Systematics
553  // Store the empirical MEC weights so we can calculate the systematics
554  // in the ShifterAndWeighter function
555  if(nu.Mode() == 1. * simb::kMEC) {
556 
557  double empiricalMECWgt = 1.;
558  double empiricalMEC_to_GENIEQE_Wgt = 1.;
559  double empiricalMEC_to_GENIERES_Wgt = 1.;
560 
561  if(tag == rwgt::kXSecCVWgt_2017){
562  // const novarwgt::EmpiricalMEC_to_GENIEQE_Wgt kEmpiricalMEC_to_GENIEQE_Wgt;
563  // const novarwgt::EmpiricalMEC_to_GENIERES_Wgt kEmpiricalMEC_to_GENIERES_Wgt;
564 
565  empiricalMECWgt = fRwgt->GetWeight(&mct, rwgt::kEmpiricalMEC_2017, &gt);
566  // empiricalMEC_to_GENIEQE_Wgt = kEmpiricalMEC_to_GENIEQE_Wgt .GetWeight(params);
567  // empiricalMEC_to_GENIERES_Wgt = kEmpiricalMEC_to_GENIERES_Wgt.GetWeight(params);
568 
569  } // end if 2017
570  else if(tag == rwgt::kXSecCVWgt_2018){
571 
572  empiricalMECWgt = fRwgt->GetWeight(&mct, rwgt::kEmpiricalMEC_2018, &gt);
573  if(vars.val_at(fnex::kTruePDG) > 0) {
574  empiricalMEC_to_GENIEQE_Wgt = fRwgt->MECq0q3ShapeSyst2018( 1, &mct, &gt);
575  empiricalMEC_to_GENIERES_Wgt = fRwgt->MECq0q3ShapeSyst2018(-1, &mct, &gt);
576  } else {
577  empiricalMEC_to_GENIEQE_Wgt = fRwgt->MECq0q3ShapeSyst2018( 1, &mct, &gt);
578  empiricalMEC_to_GENIERES_Wgt = fRwgt->MECq0q3ShapeSyst2018(-1, &mct, &gt);
579  }
580  }
581 
582  vars.set_val_at(fnex::kEmpiricalMEC_Weight, empiricalMECWgt );
583  vars.set_val_at(fnex::kEmpiricalMECtoGENIEQE_Weight, empiricalMEC_to_GENIEQE_Wgt );
584  vars.set_val_at(fnex::kEmpiricalMECtoGENIERES_Weight, empiricalMEC_to_GENIERES_Wgt);
585  //vars.set_val_at(fnex::kEmpiricalMECtoValencia_Weight, empiricalMEC_to_Valencia_Wgt);
586 
587  LOG_DEBUG("AnalysisSetupBase")
588  << "\n calculate EmpiricalMECWgt: " << empiricalMECWgt
589  << "\n calculate EmpiricalMECtoGENIEQEWgt: " << empiricalMEC_to_GENIEQE_Wgt
590  << "\n calculate EmpiricalMECtoGENIERESWgt: " << empiricalMEC_to_GENIERES_Wgt;
591  //<< "\n calculate EmpiricalMECtoValenciaWgt: " << empiricalMEC_to_Valencia_Wgt;
592  }
593  else {
597  //vars.set_val_at(fnex::kEmpiricalMECtoValencia_Weight, 1.);
598  }
599 
600  return;
601  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
double MECq0q3ShapeSyst2018(const double sigma, const simb::MCTruth *mctruth, const simb::GTruth *gtruth) const
static const unsigned char kEmpiricalMEC_Weight
Definition: VarVals.h:60
static const unsigned char kEmpiricalMECtoGENIERES_Weight
Definition: VarVals.h:62
void set_val_at(std::string const &varkey, float const &val)
Definition: VarVals.cxx:164
double GetWeight(const simb::MCTruth *mctruth, WeightType rwType, const simb::GTruth *gtruth=nullptr) const
Retrieves a weight by name. Calculates on demand, but only calculates once per neutrino.
float val_at(unsigned char const &varkey) const
Access a Var by name.
Definition: VarVals.cxx:77
static const unsigned char kEmpiricalMECtoGENIEQE_Weight
Definition: VarVals.h:61
art::ServiceHandle< rwgt::MCReweight > fRwgt
MCReweight service.
bool gt(unsigned short int a, unsigned short int b)
static const unsigned char kTruePDG
Definition: VarVals.h:43
void fnex::AnalysisSetupBase::FindRPAWeights ( fnex::MCVarVals vars,
simb::MCTruth const &  mct,
simb::GTruth const &  gt,
rwgt::WeightType const &  tag 
) const
protectedinherited

Definition at line 463 of file AnalysisSetupBase.cxx.

References std::abs(), fnex::AnalysisSetupBase::fRwgt, fnex::MCVarVals::GENIEWgts(), simb::MCTruth::GetNeutrino(), rwgt::MCReweight::GetWeight(), simb::kCC, simb::kCCQE, simb::kResCCNuBarProtonPi0Pi0, simb::kResCCNuProtonPiPlus, rwgt::kRPACCQE_2017, rwgt::kRPACCQE_SA, fnex::kRPACCQE_Weight, fnex::kRPACCQEshapeEnh, fnex::kRPACCQEshapeSupp, rwgt::kRPACCRES_2017, fnex::kRPARES_Weight, fnex::kTrueCCNC, fnex::kTrueIntType, fnex::kTruePDG, rwgt::kTuftsCC_SA, rwgt::kXSecCVWgt_2017, rwgt::kXSecCVWgt_2018, LOG_DEBUG, min(), rwgt::MCReweight::RPACCQEEnhSyst2018(), rwgt::MCReweight::RPACCQESuppSyst2018(), fnex::MCVarVals::set_GENIEWgt(), fnex::MCVarVals::set_val_at(), and fnex::MCVarVals::val_at().

Referenced by fnex::AnalysisSetupBase::FillTruthVars().

467  {
468  auto nu = mct.GetNeutrino();
469 
470  // Only set the weight from the table if the interaction type is CC QE,
471  // otherwise set it to 1 so that the calculation of the weight will return
472  // 1, ie no change in weight.
473  // The RPACCQE weights for the 2018 analysis are the same as for the 2017
474  // analysis when called from MCReweight
475  if((std::abs(vars.val_at(fnex::kTruePDG)) == 12 ||
476  std::abs(vars.val_at(fnex::kTruePDG)) == 14 ) &&
477  vars.val_at(fnex::kTrueIntType) == 1. * simb::kCCQE){
478  if(tag == rwgt::kTuftsCC_SA)
480  else if(tag == rwgt::kXSecCVWgt_2017)
482  else if(tag == rwgt::kXSecCVWgt_2018)
484  else
486  }
487  else
489 
490  // 2017 and 2018 RPA Systematics are the same
491  if(tag == rwgt::kXSecCVWgt_2017 ||
493  if(vars.val_at(fnex::kTrueCCNC) == 1. * simb::kCC &&
496 
498 
499  LOG_DEBUG("AnalysisSetupBase")
500  << "calculate RPARES("
501  << nu.QSqr()
502  << ", "
503  << double(vars.val_at(fnex::kTruePDG) < 0)
504  << ") = "
505  << fRwgt->GetWeight(&mct, rwgt::kRPACCRES_2017, &gt);
506  } // end if true cc, with appropriate resonance production
507  if(vars.val_at(fnex::kTrueIntType) == 1. * simb::kCCQE){
508  vars.set_GENIEWgt(fnex::kRPACCQEshapeEnh, "-1sigma", fRwgt->RPACCQEEnhSyst2018 (-1, &mct, &gt));
509  vars.set_GENIEWgt(fnex::kRPACCQEshapeEnh, "+1sigma", fRwgt->RPACCQEEnhSyst2018 ( 1, &mct, &gt));
510  vars.set_GENIEWgt(fnex::kRPACCQEshapeSupp, "-1sigma", fRwgt->RPACCQESuppSyst2018(-1, &mct, &gt));
511  vars.set_GENIEWgt(fnex::kRPACCQEshapeSupp, "+1sigma", fRwgt->RPACCQESuppSyst2018( 1, &mct, &gt));
512 
517 
518  LOG_DEBUG("AnalysisSetupBase")
519  << "RPACCQE CV wgt: " << fRwgt->GetWeight(&mct, rwgt::kRPACCQE_2017, &gt)
520  << "\nget RPA Enh weights (" << vars.GENIEWgts(fnex::kRPACCQEshapeEnh).at(-2)
521  << ", " << vars.GENIEWgts(fnex::kRPACCQEshapeEnh).at(-1)
522  << ", " << vars.GENIEWgts(fnex::kRPACCQEshapeEnh).at(1)
523  << ", " << vars.GENIEWgts(fnex::kRPACCQEshapeEnh).at(2)
524  << ")"
525  << "\nget RPA Supp weights (" << vars.GENIEWgts(fnex::kRPACCQEshapeSupp).at(-2)
526  << ", " << vars.GENIEWgts(fnex::kRPACCQEshapeSupp).at(-1)
527  << ", " << vars.GENIEWgts(fnex::kRPACCQEshapeSupp).at(1)
528  << ", " << vars.GENIEWgts(fnex::kRPACCQEshapeSupp).at(2)
529  << ")";
530  } // end if CCQE
531  } // end if 2017
532  else{
534 
535  vars.set_GENIEWgt(fnex::kRPACCQEshapeEnh, "-1sigma", 1.);
536  vars.set_GENIEWgt(fnex::kRPACCQEshapeEnh, "+1sigma", 1.);
537  vars.set_GENIEWgt(fnex::kRPACCQEshapeSupp, "-1sigma", 1.);
538  vars.set_GENIEWgt(fnex::kRPACCQEshapeSupp, "+1sigma", 1.);
539  }
540 
541  return;
542  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
double RPACCQEEnhSyst2018(const double sigma, const simb::MCTruth *mctruth, const simb::GTruth *gtruth) const
void set_val_at(std::string const &varkey, float const &val)
Definition: VarVals.cxx:164
static const unsigned char kRPACCQEshapeSupp
Definition: VarVals.h:75
std::map< float, float > GENIEWgts(unsigned char const &varkey)
Definition: VarVals.cxx:116
charged current quasi-elastic
Definition: MCNeutrino.h:96
float abs(float number)
Definition: d0nt_math.hpp:39
static const unsigned char kRPARES_Weight
Definition: VarVals.h:58
static const unsigned char kTrueCCNC
Definition: VarVals.h:44
double GetWeight(const simb::MCTruth *mctruth, WeightType rwType, const simb::GTruth *gtruth=nullptr) const
Retrieves a weight by name. Calculates on demand, but only calculates once per neutrino.
float val_at(unsigned char const &varkey) const
Access a Var by name.
Definition: VarVals.cxx:77
static const unsigned char kTrueIntType
Definition: VarVals.h:45
resonant charged current, nu p -> l- p pi+
Definition: MCNeutrino.h:98
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
void set_GENIEWgt(unsigned char const &varkey, std::string const &sigma, float const &wgt)
Definition: VarVals.cxx:282
static const unsigned char kRPACCQE_Weight
Definition: VarVals.h:57
art::ServiceHandle< rwgt::MCReweight > fRwgt
MCReweight service.
double RPACCQESuppSyst2018(const double sigma, const simb::MCTruth *mctruth, const simb::GTruth *gtruth) const
static const unsigned char kRPACCQEshapeEnh
Definition: VarVals.h:74
bool gt(unsigned short int a, unsigned short int b)
static const unsigned char kTruePDG
Definition: VarVals.h:43
double fnex::AnalysisSetupBase::FindXSecWeight ( simb::MCTruth const &  nu,
simb::GTruth const &  gt,
rwgt::WeightType const &  tag 
) const
protectedinherited

Definition at line 396 of file AnalysisSetupBase.cxx.

References fnex::AnalysisSetupBase::FindXSecWeightSA(), fnex::AnalysisSetupBase::fRwgt, fnex::AnalysisSetupBase::fXSecWgtScheme, rwgt::MCReweight::GetWeight(), rwgt::kTuftsCC_SA, LOG_DEBUG, and getGoodRuns4SAM::tag.

Referenced by fnex::AnalysisSetupBase::FillTruthVars().

399  {
400  LOG_DEBUG("AnalysisSetupBase")
401  << "Using weight scheme: "
402  << tag;
403 
404  if(fXSecWgtScheme == rwgt::kTuftsCC_SA ) return this->FindXSecWeightSA(mc);
405 
406  return fRwgt->GetWeight(&mc, tag, &gt);
407  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
rwgt::WeightType fXSecWgtScheme
which "Tufts" weighting to use
double GetWeight(const simb::MCTruth *mctruth, WeightType rwType, const simb::GTruth *gtruth=nullptr) const
Retrieves a weight by name. Calculates on demand, but only calculates once per neutrino.
art::ServiceHandle< rwgt::MCReweight > fRwgt
MCReweight service.
bool gt(unsigned short int a, unsigned short int b)
double FindXSecWeightSA(simb::MCTruth const &nu) const
double fnex::AnalysisSetupBase::FindXSecWeightSA ( simb::MCTruth const &  nu) const
protectedinherited

Definition at line 410 of file AnalysisSetupBase.cxx.

References wgt.

Referenced by fnex::AnalysisSetupBase::FindXSecWeight().

411  {
412  double wgt = 1.;
413 
414  // auto nu = mc.GetNeutrino();
415 
416  // // The second analysis CAF files were made without checking the EventGeneratorList
417  // // used by GENIE. For that reason, we have to explicitly set the InputVals
418  // // and call the EmpericalMECq0q3TuneWgt function here to recover their behavior.
419  // // We will actually call the DytmanMEC* weights separately too so that we can
420  // // be sure to get the correct values. The second analysis did not
421  // // check if an event was CC or NC or if it was neutrino or antineutrino.
422  // double q0 = nu.Y() * nu.Nu().E();
423  // double qmag = sqrt(nu.QSqr() + q0 * q0);
424 
425  // rwgt::DytmanMECFixXsecEdepWgt dytmanMECXsec;
426  // rwgt::DytmanMECFixItlStateWgt dytmanMECItl;
427 
428  // novarwgt::InputVals vals({{"OldGENIEMEC", true },
429  // {"qmag", qmag },
430  // {"q0", q0 },
431  // {"IsCC", true },
432  // {"IsAntiNu", nu.Nu().PdgCode() < 0},
433  // {"StruckNucleonPairPDG", nu.HitNuc() },
434  // {"Enu", nu.Nu().E() }});
435 
436  // if(nu.Mode() == simb::kMEC){
437  // // In the second analysis, an event only had to be MEC to have the
438  // // q0q3 weight applied
439  // wgt *= 0.75 * gq0q3Wgt.GetWeight(vals);
440 
441  // // In the second analysis, no checks were made on the GENIE version
442  // // to get the XsecEdep weight, here we have set the OldGENIEMEC flag
443  // // to true, which means the weight will be applied for all events
444  // wgt *= dytmanMECXsec.GetWeight(vals);
445 
446  // // In the second analysis, events were not required to be CC for the
447  // // initial state weight to be applied, but they did have to be neutrinos
448  // // in development as of 9/1/17 they were required to be CC, so we have
449  // // set the IsCC to be true for all events to replicate that behavior
450  // wgt *= dytmanMECItl.GetWeight(vals);
451  // }
452 
453  // if(nu.Mode() == simb::kDIS &&
454  // nu.W() < 1.7){
455  // rwgt::Nonres1PiWgt nonres1pi;
456  // wgt *= nonres1pi.GetWeight({{}});
457  // }
458 
459  return wgt;
460  }
const ana::Var wgt
int fnex::NuMuAnalysisSetup::FourthAnalysisQuantile ( double &  lepE,
double &  hadE,
fnex::BeamType_t const &  beamType 
) const
private

Definition at line 631 of file NuMuAnalysisSetup.cxx.

References FourthAnalysisQuantileFHC(), FourthAnalysisQuantileRHC(), fnex::kFHC, fnex::kRHC, LOG_DEBUG, and ThirdAnalysisQuantile().

Referenced by FillVars().

634  {
635 
636  if (beamType == fnex::kFHC) return FourthAnalysisQuantileFHC(lepE, hadE);
637  else if(beamType == fnex::kRHC) return FourthAnalysisQuantileRHC(lepE, hadE);
638 
639  LOG_DEBUG("NuMuAnalysisSetup")
640  << "uknown BeamType revert to Third Analysis Numu quantiles";
641  return ThirdAnalysisQuantile(lepE, hadE);
642 
643  } // fourth analysis quantiles
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
int FourthAnalysisQuantileRHC(double &lepE, double &hadE) const
int FourthAnalysisQuantileFHC(double &lepE, double &hadE) const
int ThirdAnalysisQuantile(double &recoE, double &hadEFrac) const
int fnex::NuMuAnalysisSetup::FourthAnalysisQuantileFHC ( double &  lepE,
double &  hadE 
) const
private

Definition at line 647 of file NuMuAnalysisSetup.cxx.

References LOG_DEBUG.

Referenced by FourthAnalysisQuantile().

649  {
650 
651  double recoE = lepE + hadE;
652  if( recoE == 0 ) {
653  LOG_DEBUG("NuMuAnalysisSetup")
654  << "cannot compute quantile: total energy is zero";
655  return -1;
656  }
657 
658  double hadEFrac = hadE / ( lepE + hadE );
659 
660  if( recoE >= 0 && recoE < 0.75 ){
661  if ( hadEFrac >= 0 && hadEFrac < 0.130486) return 1;
662  else if( hadEFrac >= 0.130486 && hadEFrac < 0.182243) return 2;
663  else if( hadEFrac >= 0.182243 && hadEFrac < 0.233157) return 3;
664  else if( hadEFrac >= 0.233157 && hadEFrac < 1.0) return 4;
665  }
666  else if( recoE >= 0.75 && recoE < 1 ){
667  if ( hadEFrac >= 0 && hadEFrac < 0.161791) return 1;
668  else if( hadEFrac >= 0.161791 && hadEFrac < 0.234801) return 2;
669  else if( hadEFrac >= 0.234801 && hadEFrac < 0.319627) return 3;
670  else if( hadEFrac >= 0.319627 && hadEFrac < 1.0) return 4;
671  }
672  else if( recoE >= 1 && recoE < 1.1 ){
673  if ( hadEFrac >= 0 && hadEFrac < 0.166412) return 1;
674  else if( hadEFrac >= 0.166412 && hadEFrac < 0.262069) return 2;
675  else if( hadEFrac >= 0.262069 && hadEFrac < 0.358449) return 3;
676  else if( hadEFrac >= 0.358449 && hadEFrac < 1.0) return 4;
677  }
678  else if( recoE >= 1.1 && recoE < 1.2 ){
679  if ( hadEFrac >= 0 && hadEFrac < 0.169974) return 1;
680  else if( hadEFrac >= 0.169974 && hadEFrac < 0.272577) return 2;
681  else if( hadEFrac >= 0.272577 && hadEFrac < 0.379373) return 3;
682  else if( hadEFrac >= 0.379373 && hadEFrac < 1.0) return 4;
683  }
684  else if( recoE >= 1.2 && recoE < 1.3 ){
685  if ( hadEFrac >= 0 && hadEFrac < 0.169908) return 1;
686  else if( hadEFrac >= 0.169908 && hadEFrac < 0.279256) return 2;
687  else if( hadEFrac >= 0.279256 && hadEFrac < 0.394825) return 3;
688  else if( hadEFrac >= 0.394825 && hadEFrac < 1.0) return 4;
689  }
690  else if( recoE >= 1.3 && recoE < 1.4 ){
691  if ( hadEFrac >= 0 && hadEFrac < 0.163625) return 1;
692  else if( hadEFrac >= 0.163625 && hadEFrac < 0.280182) return 2;
693  else if( hadEFrac >= 0.280182 && hadEFrac < 0.406361) return 3;
694  else if( hadEFrac >= 0.406361 && hadEFrac < 1.0) return 4;
695  }
696  else if( recoE >= 1.4 && recoE < 1.5 ){
697  if ( hadEFrac >= 0 && hadEFrac < 0.162347) return 1;
698  else if( hadEFrac >= 0.162347 && hadEFrac < 0.278829) return 2;
699  else if( hadEFrac >= 0.278829 && hadEFrac < 0.412204) return 3;
700  else if( hadEFrac >= 0.412204 && hadEFrac < 1.0) return 4;
701  }
702  else if( recoE >= 1.5 && recoE < 1.6 ){
703  if ( hadEFrac >= 0 && hadEFrac < 0.156687) return 1;
704  else if( hadEFrac >= 0.156687 && hadEFrac < 0.276816) return 2;
705  else if( hadEFrac >= 0.276816 && hadEFrac < 0.416131) return 3;
706  else if( hadEFrac >= 0.416131 && hadEFrac < 1.0) return 4;
707  }
708  else if( recoE >= 1.6 && recoE < 1.7 ){
709  if ( hadEFrac >= 0 && hadEFrac < 0.149375) return 1;
710  else if( hadEFrac >= 0.149375 && hadEFrac < 0.267368) return 2;
711  else if( hadEFrac >= 0.267368 && hadEFrac < 0.415353) return 3;
712  else if( hadEFrac >= 0.415353 && hadEFrac < 1.0) return 4;
713  }
714  else if( recoE >= 1.7 && recoE < 1.8 ){
715  if ( hadEFrac >= 0 && hadEFrac < 0.141755) return 1;
716  else if( hadEFrac >= 0.141755 && hadEFrac < 0.260274) return 2;
717  else if( hadEFrac >= 0.260274 && hadEFrac < 0.416843) return 3;
718  else if( hadEFrac >= 0.416843 && hadEFrac < 1.0) return 4;
719  }
720  else if( recoE >= 1.8 && recoE < 1.9 ){
721  if ( hadEFrac >= 0 && hadEFrac < 0.139631) return 1;
722  else if( hadEFrac >= 0.139631 && hadEFrac < 0.254985) return 2;
723  else if( hadEFrac >= 0.254985 && hadEFrac < 0.416642) return 3;
724  else if( hadEFrac >= 0.416642 && hadEFrac < 1.0) return 4;
725  }
726  else if( recoE >= 1.9 && recoE < 2 ){
727  if ( hadEFrac >= 0 && hadEFrac < 0.137587) return 1;
728  else if( hadEFrac >= 0.137587 && hadEFrac < 0.252339) return 2;
729  else if( hadEFrac >= 0.252339 && hadEFrac < 0.417286) return 3;
730  else if( hadEFrac >= 0.417286 && hadEFrac < 1.0) return 4;
731  }
732  else if( recoE >= 2 && recoE < 2.25 ){
733  if ( hadEFrac >= 0 && hadEFrac < 0.141741) return 1;
734  else if( hadEFrac >= 0.141741 && hadEFrac < 0.25798) return 2;
735  else if( hadEFrac >= 0.25798 && hadEFrac < 0.428896) return 3;
736  else if( hadEFrac >= 0.428896 && hadEFrac < 1.0) return 4;
737  }
738  else if( recoE >= 2.25 && recoE < 2.5 ){
739  if ( hadEFrac >= 0 && hadEFrac < 0.14987) return 1;
740  else if( hadEFrac >= 0.14987 && hadEFrac < 0.271534) return 2;
741  else if( hadEFrac >= 0.271534 && hadEFrac < 0.445448) return 3;
742  else if( hadEFrac >= 0.445448 && hadEFrac < 1.0) return 4;
743  }
744  else if( recoE >= 2.5 && recoE < 2.75 ){
745  if ( hadEFrac >= 0 && hadEFrac < 0.148543) return 1;
746  else if( hadEFrac >= 0.148543 && hadEFrac < 0.27805) return 2;
747  else if( hadEFrac >= 0.27805 && hadEFrac < 0.453015) return 3;
748  else if( hadEFrac >= 0.453015 && hadEFrac < 1.0) return 4;
749  }
750  else if( recoE >= 2.75 && recoE < 3 ){
751  if ( hadEFrac >= 0 && hadEFrac < 0.14472) return 1;
752  else if( hadEFrac >= 0.14472 && hadEFrac < 0.279848) return 2;
753  else if( hadEFrac >= 0.279848 && hadEFrac < 0.460109) return 3;
754  else if( hadEFrac >= 0.460109 && hadEFrac < 1.0) return 4;
755  }
756  else if( recoE >= 3 && recoE < 3.5 ){
757  if ( hadEFrac >= 0 && hadEFrac < 0.142453) return 1;
758  else if( hadEFrac >= 0.142453 && hadEFrac < 0.283427) return 2;
759  else if( hadEFrac >= 0.283427 && hadEFrac < 0.46781) return 3;
760  else if( hadEFrac >= 0.46781 && hadEFrac < 1.0) return 4;
761  }
762  else if( recoE >= 3.5 && recoE < 4 ){
763  if ( hadEFrac >= 0 && hadEFrac < 0.141182) return 1;
764  else if( hadEFrac >= 0.141182 && hadEFrac < 0.286712) return 2;
765  else if( hadEFrac >= 0.286712 && hadEFrac < 0.488556) return 3;
766  else if( hadEFrac >= 0.488556 && hadEFrac < 1.0) return 4;
767  }
768  else if( recoE >= 4 && recoE < 5 ){
769  if ( hadEFrac >= 0 && hadEFrac < 0.14111) return 1;
770  else if( hadEFrac >= 0.14111 && hadEFrac < 0.28538) return 2;
771  else if( hadEFrac >= 0.28538 && hadEFrac < 0.513618) return 3;
772  else if( hadEFrac >= 0.513618 && hadEFrac < 1.0) return 4;
773  }
774 
775  LOG_DEBUG("NuMuAnalysisSetup")
776  << "cannot compute quantile: total energy outside range";
777  return -1;
778 
779  }// fourth analysis quantiles FHC
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
int fnex::NuMuAnalysisSetup::FourthAnalysisQuantileRHC ( double &  lepE,
double &  hadE 
) const
private

Definition at line 783 of file NuMuAnalysisSetup.cxx.

References LOG_DEBUG, and LOG_VERBATIM.

Referenced by FourthAnalysisQuantile().

785  {
786 
787  double recoE = lepE + hadE;
788  if( recoE == 0 ) {
789  LOG_DEBUG("NuMuAnalysisSetup")
790  << "cannot compute quantile: total energy is zero";
791  return -1;
792  }
793 
794  double hadEFrac = hadE / ( lepE + hadE );
795 
796  LOG_VERBATIM("NuMuAnalysisSetup") << "FourthAnalysisQuantileRHC ";
797 
798  if( recoE >= 0 && recoE < 0.75 ){
799  if ( hadEFrac >= 0 && hadEFrac < 0.0980177) return 1;
800  else if( hadEFrac >= 0.0980177 && hadEFrac < 0.140308) return 2;
801  else if( hadEFrac >= 0.140308 && hadEFrac < 0.201423) return 3;
802  else if( hadEFrac >= 0.201423 && hadEFrac < 1.0) return 4;
803  }
804  else if( recoE >= 0.75 && recoE < 1 ){
805  if ( hadEFrac >= 0 && hadEFrac < 0.123805) return 1;
806  else if( hadEFrac >= 0.123805 && hadEFrac < 0.200619) return 2;
807  else if( hadEFrac >= 0.200619 && hadEFrac < 0.293348) return 3;
808  else if( hadEFrac >= 0.293348 && hadEFrac < 1.0) return 4;
809  }
810  else if( recoE >= 1 && recoE < 1.1 ){
811  if ( hadEFrac >= 0 && hadEFrac < 0.129862) return 1;
812  else if( hadEFrac >= 0.129862 && hadEFrac < 0.220222) return 2;
813  else if( hadEFrac >= 0.220222 && hadEFrac < 0.322866) return 3;
814  else if( hadEFrac >= 0.322866 && hadEFrac < 1.0) return 4;
815  }
816  else if( recoE >= 1.1 && recoE < 1.2 ){
817  if ( hadEFrac >= 0 && hadEFrac < 0.130338) return 1;
818  else if( hadEFrac >= 0.130338 && hadEFrac < 0.225444) return 2;
819  else if( hadEFrac >= 0.225444 && hadEFrac < 0.331177) return 3;
820  else if( hadEFrac >= 0.331177 && hadEFrac < 1.0) return 4;
821  }
822  else if( recoE >= 1.2 && recoE < 1.3 ){
823  if ( hadEFrac >= 0 && hadEFrac < 0.126173) return 1;
824  else if( hadEFrac >= 0.126173 && hadEFrac < 0.228904) return 2;
825  else if( hadEFrac >= 0.228904 && hadEFrac < 0.339258) return 3;
826  else if( hadEFrac >= 0.339258 && hadEFrac < 1.0) return 4;
827  }
828  else if( recoE >= 1.3 && recoE < 1.4 ){
829  if ( hadEFrac >= 0 && hadEFrac < 0.116904) return 1;
830  else if( hadEFrac >= 0.116904 && hadEFrac < 0.215843) return 2;
831  else if( hadEFrac >= 0.215843 && hadEFrac < 0.339619) return 3;
832  else if( hadEFrac >= 0.339619 && hadEFrac < 1.0) return 4;
833  }
834  else if( recoE >= 1.4 && recoE < 1.5 ){
835  if ( hadEFrac >= 0 && hadEFrac < 0.112387) return 1;
836  else if( hadEFrac >= 0.112387 && hadEFrac < 0.212124) return 2;
837  else if( hadEFrac >= 0.212124 && hadEFrac < 0.342516) return 3;
838  else if( hadEFrac >= 0.342516 && hadEFrac < 1.0) return 4;
839  }
840  else if( recoE >= 1.5 && recoE < 1.6 ){
841  if ( hadEFrac >= 0 && hadEFrac < 0.108383) return 1;
842  else if( hadEFrac >= 0.108383 && hadEFrac < 0.203952) return 2;
843  else if( hadEFrac >= 0.203952 && hadEFrac < 0.340037) return 3;
844  else if( hadEFrac >= 0.340037 && hadEFrac < 1.0) return 4;
845  }
846  else if( recoE >= 1.6 && recoE < 1.7 ){
847  if ( hadEFrac >= 0 && hadEFrac < 0.102762) return 1;
848  else if( hadEFrac >= 0.102762 && hadEFrac < 0.196337) return 2;
849  else if( hadEFrac >= 0.196337 && hadEFrac < 0.337452) return 3;
850  else if( hadEFrac >= 0.337452 && hadEFrac < 1.0) return 4;
851  }
852  else if( recoE >= 1.7 && recoE < 1.8 ){
853  if ( hadEFrac >= 0 && hadEFrac < 0.099041) return 1;
854  else if( hadEFrac >= 0.099041 && hadEFrac < 0.188072) return 2;
855  else if( hadEFrac >= 0.188072 && hadEFrac < 0.330353) return 3;
856  else if( hadEFrac >= 0.330353 && hadEFrac < 1.0) return 4;
857  }
858  else if( recoE >= 1.8 && recoE < 1.9 ){
859  if ( hadEFrac >= 0 && hadEFrac < 0.0965183) return 1;
860  else if( hadEFrac >= 0.0965183 && hadEFrac < 0.184834) return 2;
861  else if( hadEFrac >= 0.184834 && hadEFrac < 0.329758) return 3;
862  else if( hadEFrac >= 0.329758 && hadEFrac < 1.0) return 4;
863  }
864  else if( recoE >= 1.9 && recoE < 2 ){
865  if ( hadEFrac >= 0 && hadEFrac < 0.0938676) return 1;
866  else if( hadEFrac >= 0.0938676 && hadEFrac < 0.180883) return 2;
867  else if( hadEFrac >= 0.180883 && hadEFrac < 0.325615) return 3;
868  else if( hadEFrac >= 0.325615 && hadEFrac < 1.0) return 4;
869  }
870  else if( recoE >= 2 && recoE < 2.25 ){
871  if ( hadEFrac >= 0 && hadEFrac < 0.095024) return 1;
872  else if( hadEFrac >= 0.095024 && hadEFrac < 0.183673) return 2;
873  else if( hadEFrac >= 0.183673 && hadEFrac < 0.332639) return 3;
874  else if( hadEFrac >= 0.332639 && hadEFrac < 1.0) return 4;
875  }
876  else if( recoE >= 2.25 && recoE < 2.5 ){
877  if ( hadEFrac >= 0 && hadEFrac < 0.100057) return 1;
878  else if( hadEFrac >= 0.100057 && hadEFrac < 0.191958) return 2;
879  else if( hadEFrac >= 0.191958 && hadEFrac < 0.349162) return 3;
880  else if( hadEFrac >= 0.349162 && hadEFrac < 1.0) return 4;
881  }
882  else if( recoE >= 2.5 && recoE < 2.75 ){
883  if ( hadEFrac >= 0 && hadEFrac < 0.10519) return 1;
884  else if( hadEFrac >= 0.10519 && hadEFrac < 0.203698) return 2;
885  else if( hadEFrac >= 0.203698 && hadEFrac < 0.366594) return 3;
886  else if( hadEFrac >= 0.366594 && hadEFrac < 1.0) return 4;
887  }
888  else if( recoE >= 2.75 && recoE < 3 ){
889  if ( hadEFrac >= 0 && hadEFrac < 0.108997) return 1;
890  else if( hadEFrac >= 0.108997 && hadEFrac < 0.214948) return 2;
891  else if( hadEFrac >= 0.214948 && hadEFrac < 0.371937) return 3;
892  else if( hadEFrac >= 0.371937 && hadEFrac < 1.0) return 4;
893  }
894  else if( recoE >= 3 && recoE < 3.5 ){
895  if ( hadEFrac >= 0 && hadEFrac < 0.118638) return 1;
896  else if( hadEFrac >= 0.118638 && hadEFrac < 0.232254) return 2;
897  else if( hadEFrac >= 0.232254 && hadEFrac < 0.389001) return 3;
898  else if( hadEFrac >= 0.389001 && hadEFrac < 1.0) return 4;
899  }
900  else if( recoE >= 3.5 && recoE < 4 ){
901  if ( hadEFrac >= 0 && hadEFrac < 0.122369) return 1;
902  else if( hadEFrac >= 0.122369 && hadEFrac < 0.243635) return 2;
903  else if( hadEFrac >= 0.243635 && hadEFrac < 0.422904) return 3;
904  else if( hadEFrac >= 0.422904 && hadEFrac < 1.0) return 4;
905  }
906  else if( recoE >= 4 && recoE < 5 ){
907  if ( hadEFrac >= 0 && hadEFrac < 0.128739) return 1;
908  else if( hadEFrac >= 0.128739 && hadEFrac < 0.248248) return 2;
909  else if( hadEFrac >= 0.248248 && hadEFrac < 0.455024) return 3;
910  else if( hadEFrac >= 0.455024 && hadEFrac < 1.0) return 4;
911  }
912 
913  LOG_DEBUG("NuMuAnalysisSetup")
914  << "cannot compute quantile: total energy outside range";
915  return -1;
916 
917 
918  }// fourth analysis quantiles RHC
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
#define LOG_VERBATIM(category)
unsigned int fnex::NuMuAnalysisSetup::HighestPIDTrack ( std::vector< art::Ptr< rb::Track > > const &  sliceTracks,
art::InputTag const &  remidTag,
art::Event const &  e 
) const
private

Definition at line 973 of file NuMuAnalysisSetup.cxx.

References plot_validation_datamc::c, allTimeWatchdog::endl, cet::sqlite::max(), and rb::PID::Value().

976  {
977  int bestTrack = -1;
978  double highestPID = -2.0;
979  double bestTrackLength = 0.0;
980  double bestStartZ = std::numeric_limits<double>::max();
981  double bestADC = std::numeric_limits<double>::max();
982 
983  art::FindOneP<remid::ReMId> tracksToReMId(sliceTracks, e, remidTag);
984 
985  for(size_t c = 0; c < sliceTracks.size(); ++c){
986 
987  //geo::View_t view = sliceTracks[c]->View();
988  //Only use 3D tracks
989  if( sliceTracks[c]->NXCell() == 0 || sliceTracks[c]->NYCell() == 0 ){
990  continue;
991  } //End of loop over 2d tracks
992 
993  if(tracksToReMId.isValid()){
994  art::Ptr<remid::ReMId> trackReMId = tracksToReMId.at(c);
995 
996  // If the PID is lower, we can bail, nothing to update
997  if(trackReMId->Value() == highestPID)
998  { // If the new guy is shorter, we can bail
999 
1000  if (sliceTracks[c]->TotalLength() < bestTrackLength) continue;
1001 
1002  // If equal, need to fall to the next level, low start Z
1003  if (sliceTracks[c]->TotalLength() == bestTrackLength){
1004 
1005  // If start Z is higher, we can bail
1006  if (sliceTracks[c]->Start().Z() > bestStartZ) continue;
1007 
1008  // If equal, need to fall to the next level, total ADC
1009  // You might think this is stupid, but it's equivalent to flipping
1010  // a coin which always lands on the same side. This is good, this
1011  // will happen so rarely that you should just pee up a rope.
1012  if (sliceTracks[c]->Start().Z() == bestStartZ){
1013 
1014  if(sliceTracks[c]->TotalADC() > bestADC) continue;
1015 
1016  if(sliceTracks[c]->TotalADC() == bestADC)
1017  throw cet::exception("ReMId Matched Track Exception")<<std::endl
1018  << "Can't break a tie between two ReMId/tracks. "
1019  << "This should have never happend, but it did. "
1020  << "The numu group is going to have to reflect on this "
1021  << "and find a way to fix it." << std::endl;
1022  }
1023 
1024 
1025  }
1026  }
1027  highestPID = trackReMId->Value();
1028  bestTrackLength = sliceTracks[c]->TotalLength();
1029  bestStartZ = sliceTracks[c]->Start().Z();
1030  bestADC = sliceTracks[c]->TotalADC();
1031  bestTrack = c;
1032 
1033  }//End of valid ReMId object
1034  }//End of for loop over tracks
1035 
1036  if (bestTrack < 0)
1037  return 999;
1038 
1039  return (unsigned int)bestTrack;
1040 
1041  } // highestpidtrack
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
double Value() const
Definition: PID.h:22
Float_t e
Definition: plot.C:35
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
Definition: fwd.h:28
void fnex::NuMuAnalysisSetup::Reconfigure ( fhicl::ParameterSet const &  pset)
overrideprivatevirtual

The base version of this method just registers fnex::Shifters using RegisterShifters(). Derived classes can do other things in overrides (but should still call this base class method too).

Reimplemented from fnex::AnalysisSetupBase.

Definition at line 144 of file NuMuAnalysisSetup.cxx.

References fHadEFracQuantile, fNuMuEEstimator, fhicl::ParameterSet::get(), fnex::AnalysisSetupBase::Reconfigure(), and string.

Referenced by NuMuAnalysisSetup().

145  {
146  // call base class to get Cuts & shifters registered
147  this->AnalysisSetupBase::Reconfigure(pset);
148  fNuMuEEstimator = pset.get<std::string>("NuMuEEstimator", "SA");
149 
150  // get the selection type and determine the hadronic fraction qunatile
151  auto const& type = pset.get<std::string>("AnalysisType");
152 
153  fHadEFracQuantile = 0; // defalut for selections prior to TA
154 
155  if (type.compare("NuMu_Q1") == 0){
156  fHadEFracQuantile = 1;
157  }
158  else if(type.compare("NuMu_Q2") == 0){
159  fHadEFracQuantile = 2;
160  }
161  else if(type.compare("NuMu_Q3") == 0){
162  fHadEFracQuantile = 3;
163  }
164  else if(type.compare("NuMu_Q4") == 0){
165  fHadEFracQuantile = 4;
166  }
167  }
int fHadEFracQuantile
Hadronic Energy Fraction quantile.
virtual void Reconfigure(fhicl::ParameterSet const &pset)
enum BeamMode string
void fnex::NuMuAnalysisSetup::SecondAnalysisEnergy ( numue::NumuE const &  numue,
rb::Track const &  track,
double &  hadE,
double &  muonE 
) const
private

Definition at line 235 of file NuMuAnalysisSetup.cxx.

References fillBadChanDBTables::det, ds::DetectorService::DetId(), fDetService, fRH, numue::NumuE::HadCalE(), numue::NumuE::HadTrkE(), novadaq::cnv::kFARDET, novadaq::cnv::kNEARDET, numue::NumuE::NDTrkCalAct(), numue::NumuE::NDTrkCalCat(), numue::NumuE::NDTrkCalTran(), numue::NumuE::NDTrkLenAct(), numue::NumuE::NDTrkLenCat(), PandAna.reco_validation.add_data::offset, nova::dbi::RunHistory::RunNumber(), and rb::Track::TotalLength().

Referenced by FillRecoVars().

239  {
240  auto det = fDetService->DetId();
241  auto run = fRH->RunNumber();
242 
243  //New muon tuning values
244  double stitch1 = 342; // cm
245  double stitch2 = 520; // cm
246  double stitch3 = 1090; // cm
247  double offset = 0.1491; // GeV
248  double slope1 = 0.001951; // GeV/cm
249  double slope2 = 0.002022; // GeV/cm
250  double slope3 = 0.002067; // GeV/cm
251  double slope4 = 0.002196; // GeV/cm
252 
253  double trkLen = track.TotalLength();
254  double hadVisE = numue.HadTrkE() + numue.HadCalE();
255 
256  if(det == novadaq::cnv::kFARDET){
257  if (trkLen <= 0.0){
258  muonE = 0.;
259  hadE = 0.;
260  return;
261  }
262 
263  if (trkLen < stitch1){
264  muonE = slope1*trkLen + offset;
265  }
266  else if (trkLen < stitch2){
267  muonE = slope2*trkLen + ((slope1-slope2)*stitch1 + offset);
268  }
269  else if (trkLen <stitch3){
270  muonE = slope3*trkLen + ((slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
271  }
272  else{
273  muonE = slope4*trkLen + ((slope3-slope4)*stitch3 +(slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
274  }
275 
276  //New Tune Values with Error Rounding - Period 1
277  double stitchH1 = 0.0799; // GeV
278  double stitchH2 = 0.373; // GeV
279  double stitchH3 = 0.925; // GeV
280  double offsetH = 0.0419; // GeV
281  double slopeH1 = 0.816; // unitless
282  double slopeH2 = 2.04; // unitless
283  double slopeH3 = 1.68; // unitless
284  double slopeH4 = 2.45; // unitless
285 
286  if(run >= 17140 && run < 20753){
287  //New Tune Values with Error Rounding - Period 2
288  stitchH1 = 0.0694; // GeV
289  stitchH2 = 0.288; // GeV
290  stitchH3 = 0.930; // GeV
291  offsetH = 0.0392; // GeV
292  slopeH1 = 1.45; // unitless
293  slopeH2 = 2.13; // unitless
294  slopeH3 = 1.76; // unitless
295  slopeH4 = 2.17; // unitless
296  }
297  else if(run >= 20753){
298  //New Tune Values with Error Rounding - Period 3
299  stitchH1 = 0.0702; // GeV
300  stitchH2 = 0.321; // GeV
301  stitchH3 = 1.06; // GeV
302  offsetH = 0.0365; // GeV
303  slopeH1 = 1.34; // unitless
304  slopeH2 = 2.05; // unitless
305  slopeH3 = 1.68; // unitless
306  slopeH4 = 2.28; // unitless
307  }
308 
309  if (hadVisE < 0.0){
310  hadE = offsetH;
311  }
312  else if (hadVisE < stitchH1){
313  hadE = slopeH1*hadVisE + offsetH;
314  }
315  else if (hadVisE < stitchH2){
316  hadE = slopeH2*hadVisE + ((slopeH1-slopeH2)*stitchH1 + offsetH);
317  }
318  else if (hadVisE <stitchH3){
319  hadE = slopeH3*hadVisE + ((slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
320  }
321  else{
322  hadE = slopeH4*hadVisE + ((slopeH3-slopeH4)*stitchH3 +(slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
323  }
324  }//FD energy
325 
327  double trkLenAct = numue.NDTrkLenAct();
328  double trkLenCat = numue.NDTrkLenCat();
329 
330  if(trkLenAct <= 0.0 && trkLenCat <= 0.0){
331  muonE = 0.;
332  hadE = 0.;
333  return;
334  }
335 
336  bool actOnly = ((numue.NDTrkCalAct() > 0.0 ||
337  numue.NDTrkCalTran() > 0.0 ) &&
338  numue.NDTrkCalCat() == 0.0 );
339  bool actAndCat = ((numue.NDTrkCalAct() > 0.0 ||
340  numue.NDTrkCalTran() > 0.0 ) &&
341  numue.NDTrkCalCat() > 0.0 );
342 
343  if( !actOnly && !actAndCat){
344  muonE = 0.;
345  hadE = 0.;
346  return;
347  }
348  if( actOnly && actAndCat){
349  muonE = 0.;
350  hadE = 0.;
351  }
352 
353  if(actOnly){
354  trkLen = trkLenAct;
355  }
356 
357  if(actAndCat){
358 
359  double muonEMC = 0.0;
360 
361  //New tuning values rounded for errors
362  double offsetMC = 0.139; // GeV
363  double slope1MC = 0.00537; // GeV/cm
364 
365  muonEMC = slope1MC*trkLenCat + offsetMC;
366 
367  double effActTrkLen = 0.0;
368 
369  if(muonEMC <= 0.0){
370  muonE = 0.;
371  hadE = 0.;
372  return;
373  }
374 
375  if(muonEMC <= offset) {
376  effActTrkLen = 0.0;
377  }
378  else if(muonEMC < (slope1*stitch1 + offset)){
379  effActTrkLen = (muonEMC - offset)/slope1;
380  }
381  else if(muonEMC < (slope2*stitch2 + ((slope1-slope2)*stitch1 + offset))){
382  effActTrkLen = (muonEMC - ((slope1-slope2)*stitch1 + offset))/slope2;
383  }
384  else if(muonEMC < (slope3*stitch3 + ((slope2-slope3)*stitch2 + (slope1-slope2)*stitch1 + offset))){
385  effActTrkLen = (muonEMC - ((slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset))/slope3;
386  }
387  else{
388  effActTrkLen = (muonEMC - ((slope3-slope4)*stitch3 + (slope2-slope3)*stitch2 + (slope1-slope2)*stitch1 + offset))/slope4;
389  }
390 
391  trkLen = trkLenAct + effActTrkLen;
392 
393  }//End of actAndCat
394 
395  //Now get final muon energy from total "active" track length
396  if(trkLen < stitch1){
397  muonE = slope1*trkLen + offset;
398  }
399  else if (trkLen < stitch2){
400  muonE = slope2*trkLen + ((slope1-slope2)*stitch1 + offset);
401  }
402  else if (trkLen <stitch3){
403  muonE = slope3*trkLen + ((slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
404  }
405  else{
406  muonE = slope4*trkLen + ((slope3-slope4)*stitch3 +(slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
407  }
408 
409 
410  //New Tune Values with Error Rounding
411  double stitchH1 = 0.0394; // GeV
412  double stitchH2 = 0.362; // GeV
413  double stitchH3 = 0.955; // GeV
414  double offsetH = 0.0534; // GeV
415  double slopeH1 = 0.956; // unitless
416  double slopeH2 = 2.14; // unitless
417  double slopeH3 = 1.78; // unitless
418  double slopeH4 = 2.00; // unitless
419 
420  if (hadVisE < 0.0){
421  hadE = offsetH;
422  }
423  else if (hadVisE < stitchH1){
424  hadE = slopeH1*hadVisE + offsetH;
425  }
426  else if (hadVisE < stitchH2){
427  hadE = slopeH2*hadVisE + ((slopeH1-slopeH2)*stitchH1 + offsetH);
428  }
429  else if (hadVisE <stitchH3){
430  hadE = slopeH3*hadVisE + ((slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
431  }
432  else{
433  hadE = slopeH4*hadVisE + ((slopeH3-slopeH4)*stitchH3 +(slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
434  }
435 
436  }//ND energy
437 
438  // START TERRIBLE HACK FOR SYSTEMATICS
439 // if(tMuE > 0) muonE = tMuE;
440 // if(tHadE > 0) hadE = tHadE;
441 // double nuE = muonE + hadE;
442 // if(tNuEScale > 0) nuE *= tNuEScale;
443 // return nuE;
444 
445  }
novadaq::cnv::DetId DetId() const
What detector are we in?
Energy estimators for CC events.
Definition: FillEnergies.h:7
Definition: event.h:19
Far Detector at Ash River, MN.
art::ServiceHandle< nova::dbi::RunHistoryService > fRH
RunHistory service.
Near Detector in the NuMI cavern.
Definition: run.py:1
art::ServiceHandle< ds::DetectorService > fDetService
which detector?
int RunNumber() const
Definition: RunHistory.h:374
void fnex::NuMuAnalysisSetup::ThirdAnalysisEnergy ( numue::NumuE const &  numue,
rb::Track const &  track,
double &  hadE,
double &  muonE,
bool  isMC,
fnex::BeamType_t const &  beamType 
) const
private

Definition at line 448 of file NuMuAnalysisSetup.cxx.

References fillBadChanDBTables::det, ds::DetectorService::DetId(), fDetService, fRH, numue::NumuE::HadCalE(), numue::NumuE::HadTrkE(), make_mec_shifts_plots::isRHC, novadaq::cnv::kFARDET, novadaq::cnv::kNEARDET, fnex::kRHC, LOG_DEBUG, numue::NumuE::NDTrkLenAct(), numue::NumuE::NDTrkLenCat(), predict_prod3_fd_had_energy(), predict_prod3_fd_muon_energy(), predict_prod3_nd_act_energy(), predict_prod3_nd_cat_energy(), predict_prod3_nd_had_energy(), nova::dbi::RunHistory::RunNumber(), and rb::Track::TotalLength().

Referenced by FillRecoVars().

454  {
455  auto det = fDetService->DetId();
456  auto run = fRH->RunNumber();
457 
458  double hadVisE = numue.HadTrkE() + numue.HadCalE();
459 
460  bool isRHC = false; // FHC mode
461  if(beamType == fnex::kRHC) isRHC = true ; // RHC mode
462 
463 
465  {
466 
467  double trkLen = track.TotalLength();
468 
469  // if(isMC && !isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "MC, not RHC";
470  // if(!isMC && !isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "Not MC, not RHC";
471 
472  muonE = predict_prod3_fd_muon_energy(trkLen, run, isMC, isRHC);
473  hadE = predict_prod3_fd_had_energy(hadVisE, run, isMC, isRHC);
474 
475  LOG_DEBUG("NuMuAnalysisSetup") << "LepE = " << muonE << "\t HadE = " << hadE;
476  }
477 
479  {
480 
481  double trkLenAct = numue.NDTrkLenAct();
482  double trkLenCat = numue.NDTrkLenCat();
483 
484  // if(isMC && !isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "MC, not RHC";
485  // if(!isMC && !isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "Not MC, not RHC";
486 
487  muonE = predict_prod3_nd_act_energy(trkLenAct, isRHC) + predict_prod3_nd_cat_energy(trkLenCat, isRHC);
488  hadE = predict_prod3_nd_had_energy(hadVisE, isRHC);
489  }
490 
491  }// end TA energy estimator
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
novadaq::cnv::DetId DetId() const
What detector are we in?
Energy estimators for CC events.
Definition: FillEnergies.h:7
Definition: event.h:19
double predict_prod3_fd_muon_energy(double trkLen, int run, bool ismc, bool isRHC)
Far Detector at Ash River, MN.
art::ServiceHandle< nova::dbi::RunHistoryService > fRH
RunHistory service.
Near Detector in the NuMI cavern.
Definition: run.py:1
art::ServiceHandle< ds::DetectorService > fDetService
which detector?
double predict_prod3_nd_cat_energy(double ndtrklencat, bool isRHC)
double predict_prod3_nd_had_energy(double hadVisE, bool isRHC)
int RunNumber() const
Definition: RunHistory.h:374
double predict_prod3_fd_had_energy(double hadVisE, int run, bool ismc, bool isRHC)
double predict_prod3_nd_act_energy(double ndtrklenact, bool isRHC)
int fnex::NuMuAnalysisSetup::ThirdAnalysisQuantile ( double &  recoE,
double &  hadEFrac 
) const
private

Definition at line 495 of file NuMuAnalysisSetup.cxx.

References LOG_DEBUG.

Referenced by FourthAnalysisQuantile().

497  {
498 
499  double recoE = lepE + hadE;
500  if( recoE == 0 ) {
501  LOG_DEBUG("NuMuAnalysisSetup")
502  << "cannot compute quantile: total energy is zero";
503  return -1;
504  }
505 
506  double hadEFrac = hadE / ( lepE + hadE );
507 
508  if( recoE >= 0 && recoE < 0.75 ){
509  if ( hadEFrac >= 0 && hadEFrac < 0.160864) return 1;
510  else if( hadEFrac >= 0.160864 && hadEFrac < 0.231946) return 2;
511  else if( hadEFrac >= 0.231946 && hadEFrac < 0.294404) return 3;
512  else if( hadEFrac >= 0.294404 && hadEFrac <= 1 ) return 4;
513  }
514  else if( recoE >= 0.75 && recoE < 1 ){
515  if ( hadEFrac >= 0 && hadEFrac < 0.201787) return 1;
516  else if( hadEFrac >= 0.201787 && hadEFrac < 0.290059) return 2;
517  else if( hadEFrac >= 0.290059 && hadEFrac < 0.376756) return 3;
518  else if( hadEFrac >= 0.376756 && hadEFrac <= 1 ) return 4;
519  }
520  else if( recoE >= 1 && recoE < 1.1 ){
521  if ( hadEFrac >= 0 && hadEFrac < 0.209212) return 1;
522  else if( hadEFrac >= 0.209212 && hadEFrac < 0.306823) return 2;
523  else if( hadEFrac >= 0.306823 && hadEFrac < 0.404893) return 3;
524  else if( hadEFrac >= 0.404893 && hadEFrac <= 1 ) return 4;
525  }
526  else if( recoE >= 1.1 && recoE < 1.2 ){
527  if ( hadEFrac >= 0 && hadEFrac < 0.214098) return 1;
528  else if( hadEFrac >= 0.214098 && hadEFrac < 0.314586) return 2;
529  else if( hadEFrac >= 0.314586 && hadEFrac < 0.417883) return 3;
530  else if( hadEFrac >= 0.417883 && hadEFrac <= 1 ) return 4;
531  }
532  else if( recoE >= 1.2 && recoE < 1.3 ){
533  if ( hadEFrac >= 0 && hadEFrac < 0.209203) return 1;
534  else if( hadEFrac >= 0.209203 && hadEFrac < 0.316993) return 2;
535  else if( hadEFrac >= 0.316993 && hadEFrac < 0.428731) return 3;
536  else if( hadEFrac >= 0.428731 && hadEFrac <= 1 ) return 4;
537  }
538  else if( recoE >= 1.3 && recoE < 1.4 ){
539  if ( hadEFrac >= 0 && hadEFrac < 0.203099) return 1;
540  else if( hadEFrac >= 0.203099 && hadEFrac < 0.315997) return 2;
541  else if( hadEFrac >= 0.315997 && hadEFrac < 0.432587) return 3;
542  else if( hadEFrac >= 0.432587 && hadEFrac <= 1 ) return 4;
543  }
544  else if( recoE >= 1.4 && recoE < 1.5 ){
545  if ( hadEFrac >= 0 && hadEFrac < 0.198472) return 1;
546  else if( hadEFrac >= 0.198472 && hadEFrac < 0.31283 ) return 2;
547  else if( hadEFrac >= 0.31283 && hadEFrac < 0.442388) return 3;
548  else if( hadEFrac >= 0.442388 && hadEFrac <= 1 ) return 4;
549  }
550  else if( recoE >= 1.5 && recoE < 1.6 ){
551  if ( hadEFrac >= 0 && hadEFrac < 0.188711) return 1;
552  else if( hadEFrac >= 0.188711 && hadEFrac < 0.302839) return 2;
553  else if( hadEFrac >= 0.302839 && hadEFrac < 0.442194) return 3;
554  else if( hadEFrac >= 0.442194 && hadEFrac <= 1 ) return 4;
555  }
556  else if( recoE >= 1.6 && recoE < 1.7 ){
557  if ( hadEFrac >= 0 && hadEFrac < 0.180827) return 1;
558  else if( hadEFrac >= 0.180827 && hadEFrac < 0.293246) return 2;
559  else if( hadEFrac >= 0.293246 && hadEFrac < 0.443012) return 3;
560  else if( hadEFrac >= 0.443012 && hadEFrac <= 1 ) return 4;
561  }
562  else if( recoE >= 1.7 && recoE < 1.8 ){
563  if ( hadEFrac >= 0 && hadEFrac < 0.172282) return 1;
564  else if( hadEFrac >= 0.172282 && hadEFrac < 0.283574) return 2;
565  else if( hadEFrac >= 0.283574 && hadEFrac < 0.441167) return 3;
566  else if( hadEFrac >= 0.441167 && hadEFrac <= 1 ) return 4;
567  }
568  else if( recoE >= 1.8 && recoE < 1.9 ){
569  if ( hadEFrac >= 0 && hadEFrac < 0.16746) return 1;
570  else if( hadEFrac >= 0.16746 && hadEFrac < 0.27588) return 2;
571  else if( hadEFrac >= 0.27588 && hadEFrac < 0.43916) return 3;
572  else if( hadEFrac >= 0.43916 && hadEFrac <= 1 ) return 4;
573  }
574  else if( recoE >= 1.9 && recoE < 2 ){
575  if ( hadEFrac >= 0 && hadEFrac < 0.1646 ) return 1;
576  else if( hadEFrac >= 0.1646 && hadEFrac < 0.272764) return 2;
577  else if( hadEFrac >= 0.272764 && hadEFrac < 0.444039) return 3;
578  else if( hadEFrac >= 0.444039 && hadEFrac <= 1 ) return 4;
579  }
580  else if( recoE >= 2 && recoE < 2.25 ){
581  if ( hadEFrac >= 0 && hadEFrac < 0.165357) return 1;
582  else if( hadEFrac >= 0.165357 && hadEFrac < 0.27372) return 2;
583  else if( hadEFrac >= 0.27372 && hadEFrac < 0.451054) return 3;
584  else if( hadEFrac >= 0.451054 && hadEFrac <= 1 ) return 4;
585  }
586  else if( recoE >= 2.25 && recoE < 2.5 ){
587  if ( hadEFrac >= 0 && hadEFrac < 0.169381) return 1;
588  else if( hadEFrac >= 0.169381 && hadEFrac < 0.282544) return 2;
589  else if( hadEFrac >= 0.282544 && hadEFrac < 0.460158) return 3;
590  else if( hadEFrac >= 0.460158 && hadEFrac <= 1 ) return 4;
591  }
592  else if( recoE >= 2.5 && recoE < 2.75 ){
593  if ( hadEFrac >= 0 && hadEFrac < 0.165923) return 1;
594  else if( hadEFrac >= 0.165923 && hadEFrac < 0.283755) return 2;
595  else if( hadEFrac >= 0.283755 && hadEFrac < 0.456392) return 3;
596  else if( hadEFrac >= 0.456392 && hadEFrac <= 1 ) return 4;
597  }
598  else if( recoE >= 2.75 && recoE < 3 ){
599  if ( hadEFrac >= 0 && hadEFrac < 0.158863) return 1;
600  else if( hadEFrac >= 0.158863 && hadEFrac < 0.283625) return 2;
601  else if( hadEFrac >= 0.283625 && hadEFrac < 0.45765 ) return 3;
602  else if( hadEFrac >= 0.45765 && hadEFrac <= 1 ) return 4;
603  }
604  else if( recoE >= 3 && recoE < 3.5 ){
605  if ( hadEFrac >= 0 && hadEFrac < 0.157646) return 1;
606  else if( hadEFrac >= 0.157646 && hadEFrac < 0.290466) return 2;
607  else if( hadEFrac >= 0.290466 && hadEFrac < 0.468987) return 3;
608  else if( hadEFrac >= 0.468987 && hadEFrac <= 1 ) return 4;
609  }
610  else if( recoE >= 3.5 && recoE < 4 ){
611  if ( hadEFrac >= 0 && hadEFrac < 0.153759) return 1;
612  else if( hadEFrac >= 0.153759 && hadEFrac < 0.295047) return 2;
613  else if( hadEFrac >= 0.295047 && hadEFrac < 0.495104) return 3;
614  else if( hadEFrac >= 0.495104 && hadEFrac <= 1 ) return 4;
615  }
616  else if( recoE >= 4 && recoE < 5 ){
617  if ( hadEFrac >= 0 && hadEFrac < 0.154679) return 1;
618  else if( hadEFrac >= 0.154679 && hadEFrac < 0.299486) return 2;
619  else if( hadEFrac >= 0.299486 && hadEFrac < 0.542474) return 3;
620  else if( hadEFrac >= 0.542474 && hadEFrac <= 1 ) return 4;
621  }
622 
623  LOG_DEBUG("NuMuAnalysisSetup")
624  << "cannot compute quantile: total energy outside range";
625  return -1;
626 
627  }// quantiles
#define LOG_DEBUG(stream)
Definition: Messenger.h:149

Member Data Documentation

art::ServiceHandle<ds::DetectorService> fnex::NuMuAnalysisSetup::fDetService
private

which detector?

Definition at line 101 of file NuMuAnalysisSetup.h.

Referenced by Analysis2018Energy(), FillRecoVars(), SecondAnalysisEnergy(), and ThirdAnalysisEnergy().

int fnex::NuMuAnalysisSetup::fHadEFracQuantile
private

Hadronic Energy Fraction quantile.

Definition at line 98 of file NuMuAnalysisSetup.h.

Referenced by FillVars(), and Reconfigure().

std::string fnex::NuMuAnalysisSetup::fNuMuEEstimator
private

Definition at line 54 of file NuMuAnalysisSetup.h.

Referenced by FillRecoVars(), and Reconfigure().

art::ServiceHandle<nova::dbi::RunHistoryService> fnex::NuMuAnalysisSetup::fRH
private

RunHistory service.

Definition at line 100 of file NuMuAnalysisSetup.h.

Referenced by Analysis2018Energy(), SecondAnalysisEnergy(), and ThirdAnalysisEnergy().

art::ServiceHandle<rwgt::MCReweight> fnex::AnalysisSetupBase::fRwgt
protectedinherited
art::ServiceHandle<art::TFileService> fnex::AnalysisSetupBase::fTFS
protectedinherited

TFileService.

Definition at line 110 of file AnalysisSetupBase.h.

rwgt::WeightType fnex::AnalysisSetupBase::fXSecWgtScheme
protectedinherited

The documentation for this class was generated from the following files: