Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
fnex::AnalysisSetupBase Class Referenceabstract

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-11-28/FNEX/core/AnalysisSetupBase.h"

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

Public Member Functions

virtual ~AnalysisSetupBase ()=default
 constructs a VarList with Vars filled from each slice in this Event More...
 
virtual 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 =0
 
virtual void Reconfigure (fhicl::ParameterSet const &pset)
 

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...
 

Detailed Description

Definition at line 59 of file AnalysisSetupBase.h.

Constructor & Destructor Documentation

virtual fnex::AnalysisSetupBase::~AnalysisSetupBase ( )
virtualdefault

constructs a VarList with Vars filled from each slice in this Event

Member Function Documentation

void fnex::AnalysisSetupBase::FillGENIEWeightVars ( fnex::MCVarVals vars,
rwgt::GENIEReweightTable const &  rwgtTable 
) const
protected

Definition at line 47 of file AnalysisSetupBase.cxx.

References 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 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
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
protected

Definition at line 211 of file AnalysisSetupBase.cxx.

References evt, FillGENIEWeightVars(), FindDISWeight(), FindMECWeight(), FindRPAWeights(), 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, 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 fnex::NuMuAnalysisSetup::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
virtual void fnex::AnalysisSetupBase::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
pure virtual
void fnex::AnalysisSetupBase::FindDISWeight ( fnex::MCVarVals vars,
simb::MCTruth const &  mct,
simb::GTruth const &  gt,
rwgt::WeightType const &  tag 
) const
protected

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 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
protected

Definition at line 545 of file AnalysisSetupBase.cxx.

References 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 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
protected

Definition at line 463 of file AnalysisSetupBase.cxx.

References std::abs(), 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 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
protected

Definition at line 396 of file AnalysisSetupBase.cxx.

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

Referenced by 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
protected

Definition at line 410 of file AnalysisSetupBase.cxx.

References wgt.

Referenced by 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
void fnex::AnalysisSetupBase::Reconfigure ( fhicl::ParameterSet const &  pset)
virtual

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 in fnex::NuEAnalysisSetup, and fnex::NuMuAnalysisSetup.

Definition at line 35 of file AnalysisSetupBase.cxx.

References fXSecWgtScheme, fhicl::ParameterSet::get(), rwgt::kTuftsCC_SA, rwgt::kXSecCVWgt_2017, rwgt::kXSecCVWgt_2018, and string.

Referenced by fnex::NuMuAnalysisSetup::Reconfigure(), and fnex::NuEAnalysisSetup::Reconfigure().

36  {
37  auto wgtToUse = pset.get<std::string>("XSecWgtScheme", "SA");
38 
39  if(wgtToUse.compare("SA") == 0) fXSecWgtScheme = rwgt::kTuftsCC_SA;
40  else if(wgtToUse.compare("2017") == 0) fXSecWgtScheme = rwgt::kXSecCVWgt_2017;
41  else if(wgtToUse.compare("2018") == 0) fXSecWgtScheme = rwgt::kXSecCVWgt_2018;
42 
43  return;
44  }
rwgt::WeightType fXSecWgtScheme
which "Tufts" weighting to use
enum BeamMode string

Member Data Documentation

art::ServiceHandle<rwgt::MCReweight> fnex::AnalysisSetupBase::fRwgt
protected

MCReweight service.

Definition at line 111 of file AnalysisSetupBase.h.

Referenced by FindMECWeight(), FindRPAWeights(), and FindXSecWeight().

art::ServiceHandle<art::TFileService> fnex::AnalysisSetupBase::fTFS
protected

TFileService.

Definition at line 110 of file AnalysisSetupBase.h.

rwgt::WeightType fnex::AnalysisSetupBase::fXSecWgtScheme
protected

which "Tufts" weighting to use

Definition at line 112 of file AnalysisSetupBase.h.

Referenced by FillGENIEWeightVars(), FillTruthVars(), FindXSecWeight(), and Reconfigure().


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