Public Member Functions | Public Attributes | Private Attributes | List of all members
NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter Class Reference

Reweight a chain of interactions that are covered by the NuMI target K/pi ratios measured by MIPP. More...

#include <MIPPNumiKaonYieldsReweighter.h>

Inheritance diagram for NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter:
NeutrinoFluxReweight::IInteractionChainReweighting

Public Member Functions

 MIPPNumiKaonYieldsReweighter (int iuniv, const ParameterTable &cv_pars, const ParameterTable &univ_pars)
 
virtual ~MIPPNumiKaonYieldsReweighter ()
 
virtual std::vector< bool > canReweight (const InteractionChainData &aa)
 Look through the InteractionChainData input and identify those Interactions that can be reweighted as part of a chain. We return a vector indicating which elements will be assigned a weight by calculateWeight. More...
 
virtual double calculateWeight (const InteractionChainData &aa)
 calculate a weight for this interaction chain given the central value parameters and the parameters for this universe. The weight is something like: f(cv)/f(MC) * f(univ)/f(cv) where cv in this case corresponds to the best value of the parameter, given the data. If univ_pars=cv_pars then we are calculating a central value weight. Note, canReweight() should be called to determine which elements of the chain are covered by the weight returned by calculateWeight() More...
 
- Public Member Functions inherited from NeutrinoFluxReweight::IInteractionChainReweighting
virtual ~IInteractionChainReweighting ()
 

Public Attributes

const ParameterTable & cvPars
 
const ParameterTable & univPars
 

Private Attributes

int iUniv
 
float prt_no_inter
 
std::vector< float > vbin_datacv_pip
 
std::vector< float > vbin_datacv_pim
 
std::vector< float > vbin_datasys_pip
 
std::vector< float > vbin_datasys_pim
 
std::vector< float > vbin_datasta_pip
 
std::vector< float > vbin_datasta_pim
 
std::vector< float > vbin_datacv_kap_pip
 
std::vector< float > vbin_datacv_kam_pim
 
std::vector< float > vbin_datasys_kap_pip
 
std::vector< float > vbin_datasys_kam_pim
 
std::vector< float > vbin_datasta_kap_pip
 
std::vector< float > vbin_datasta_kam_pim
 
float aux_par
 

Detailed Description

Reweight a chain of interactions that are covered by the NuMI target K/pi ratios measured by MIPP.

Definition at line 15 of file MIPPNumiKaonYieldsReweighter.h.

Constructor & Destructor Documentation

NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::MIPPNumiKaonYieldsReweighter ( int  iuniv,
const ParameterTable &  cv_pars,
const ParameterTable &  univ_pars 
)

The constructor. Note, we pass central value and single universe parameters in this constructor only. There is thus a 1 to 1 correspondence between an instance of this class and a given universe.

Definition at line 11 of file MIPPNumiKaonYieldsReweighter.cpp.

11  : cvPars(cv_pars), univPars(univ_pars), iUniv(iuniv){
12 
13  MIPPNumiYieldsBins* MIPPbins = MIPPNumiYieldsBins::getInstance();
14  vbin_datacv_pip.reserve(MIPPbins->GetNbins_pip_MIPPNuMI());
15  vbin_datacv_pim.reserve(MIPPbins->GetNbins_pim_MIPPNuMI());
16  vbin_datasys_pip.reserve(MIPPbins->GetNbins_pip_MIPPNuMI());
17  vbin_datasys_pim.reserve(MIPPbins->GetNbins_pim_MIPPNuMI());
18  vbin_datasta_pip.reserve(MIPPbins->GetNbins_pip_MIPPNuMI());
19  vbin_datasta_pim.reserve(MIPPbins->GetNbins_pim_MIPPNuMI());
20 
21  vbin_datacv_kap_pip.reserve(MIPPbins->GetNbins_K_MIPPNuMI());
22  vbin_datacv_kam_pim.reserve(MIPPbins->GetNbins_K_MIPPNuMI());
23  vbin_datasys_kap_pip.reserve(MIPPbins->GetNbins_K_MIPPNuMI());
24  vbin_datasys_kam_pim.reserve(MIPPbins->GetNbins_K_MIPPNuMI());
25  vbin_datasta_kap_pip.reserve(MIPPbins->GetNbins_K_MIPPNuMI());
26  vbin_datasta_kam_pim.reserve(MIPPbins->GetNbins_K_MIPPNuMI());
27 
28  //const boost::interprocess::flat_map<std::string, double>& univ_table = univPars.getMap();
29  //const boost::interprocess::flat_map<std::string, double>& cv_table = cvPars.getMap();
30  prt_no_inter = univPars.getParameterValue("prt_no_interacting");
31 
32  char namepar[100];
33  for(int ii=0;ii<MIPPbins->GetNbins_pip_MIPPNuMI();ii++){
34  sprintf(namepar,"MIPP_NuMI_%s_sys_%d","pip",ii);
35  double data_cv = cvPars.getParameterValue(std::string(namepar));
36  double data_sys = univPars.getParameterValue(std::string(namepar));
37  sprintf(namepar,"MIPP_NuMI_%s_stats_%d","pip",ii);
38  double data_sta = univPars.getParameterValue(std::string(namepar));
39  data_sys /= (1.0-prt_no_inter);
40  data_sta /= (1.0-prt_no_inter);
41  data_cv /= (1.0-prt_no_inter);
42  vbin_datacv_pip.push_back(data_cv);
43  vbin_datasys_pip.push_back(data_sys);
44  vbin_datasta_pip.push_back(data_sta);
45  }
46 
47  for(int ii=0;ii<MIPPbins->GetNbins_pim_MIPPNuMI();ii++){
48  sprintf(namepar,"MIPP_NuMI_%s_sys_%d","pim",ii);
49  double data_cv = cvPars.getParameterValue(std::string(namepar));
50  double data_sys = univPars.getParameterValue(std::string(namepar));
51  sprintf(namepar,"MIPP_NuMI_%s_stats_%d","pim",ii);
52  double data_sta = univPars.getParameterValue(std::string(namepar));
53  data_sys /= (1.0-prt_no_inter);
54  data_sta /= (1.0-prt_no_inter);
55  data_cv /= (1.0-prt_no_inter);
56  vbin_datacv_pim.push_back(data_cv);
57  vbin_datasys_pim.push_back(data_sys);
58  vbin_datasta_pim.push_back(data_sta);
59  }
60  for(int ii=0;ii<MIPPbins->GetNbins_K_MIPPNuMI();ii++){
61  sprintf(namepar,"MIPP_NuMI_%s_sys_%d","kap_pip",ii);
62  double data_cv = cvPars.getParameterValue(std::string(namepar));
63  double data_sys = univPars.getParameterValue(std::string(namepar));
64  sprintf(namepar,"MIPP_NuMI_%s_stats_%d","kap_pip",ii);
65  double data_sta = univPars.getParameterValue(std::string(namepar));
66  vbin_datacv_kap_pip.push_back(data_cv);
67  vbin_datasys_kap_pip.push_back(data_sys);
68  vbin_datasta_kap_pip.push_back(data_sta);
69 
70  sprintf(namepar,"MIPP_NuMI_%s_sys_%d","kam_pim",ii);
71  data_cv = cvPars.getParameterValue(std::string(namepar));
72  data_sys = univPars.getParameterValue(std::string(namepar));
73  sprintf(namepar,"MIPP_NuMI_%s_stats_%d","kam_pim",ii);
74  data_sta = univPars.getParameterValue(std::string(namepar));
75  vbin_datacv_kam_pim.push_back(data_cv);
76  vbin_datasys_kam_pim.push_back(data_sys);
77  vbin_datasta_kam_pim.push_back(data_sta);
78  }
79  aux_par = univPars.getParameterValue("aux_parameter");
80  if(aux_par<1.e-15)aux_par = 1.0;
81 
82  }
string
std::string string
Definition: nybbler.cc:12
NeutrinoFluxReweight::ParameterTable::getParameterValue
double getParameterValue(const std::string &name) const
get the value of a parameter. throw an exception of a well defined type if we don&#39;t have it ...
Definition: ParameterTable.cpp:33
e
const double e
Definition: gUpMuFluxGen.cxx:165
NeutrinoFluxReweight::MIPPNumiYieldsBins::getInstance
static MIPPNumiYieldsBins * getInstance()
Definition: MIPPNumiYieldsBins.cpp:153
NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::~MIPPNumiKaonYieldsReweighter ( )
virtual

Definition at line 83 of file MIPPNumiKaonYieldsReweighter.cpp.

83  {
84 
85  }

Member Function Documentation

double NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::calculateWeight ( const InteractionChainData &  aa)
virtual

calculate a weight for this interaction chain given the central value parameters and the parameters for this universe. The weight is something like: f(cv)/f(MC) * f(univ)/f(cv) where cv in this case corresponds to the best value of the parameter, given the data. If univ_pars=cv_pars then we are calculating a central value weight. Note, canReweight() should be called to determine which elements of the chain are covered by the weight returned by calculateWeight()

Implements NeutrinoFluxReweight::IInteractionChainReweighting.

Definition at line 135 of file MIPPNumiKaonYieldsReweighter.cpp.

135  {
136 
137  double wgt = 1.0;
138 
139  MIPPNumiYieldsBins* MIPPbins = MIPPNumiYieldsBins::getInstance();
140  MIPPNumiMC* MCval = MIPPNumiMC::getInstance();
141  double low_value = 1.e-18;
142 
143  TargetData tar = aa.tar_info;
144 
145  //fast check:
146  if(tar.Tar_pdg != 321 && tar.Tar_pdg != -321 && tar.Tar_pdg != 130 && tar.Tar_pdg != 310)return aux_par;
147  if(tar.Pz<20.0 || tar.Pz>80.0 || tar.Pt>2.0)return aux_par;
148  int binID = MIPPbins->BinID(tar.Pz,tar.Pt,tar.Tar_pdg);
149  if(binID<0){
150  return aux_par;
151  }
152  //Now looking for pions bin ID:
153  int pip_bin = MIPPbins->BinID(tar.Pz,tar.Pt, 211);
154  int pim_bin = MIPPbins->BinID(tar.Pz,tar.Pt,-211);
155  if(tar.Tar_pdg == 321 && pip_bin<0)return aux_par;
156  if(tar.Tar_pdg ==-321 && pim_bin<0)return aux_par;
157  if(tar.Tar_pdg == 130 || tar.Tar_pdg == 310){
158  if(pip_bin<0 || pim_bin<0)return aux_par;
159  }
160 
161  //Now, looking for the MC value:
162  double binC = MCval->getMCval(tar.Pz,tar.Pt,tar.Tar_pdg);
163  if(binC<low_value){
164  //std::cout<<"LOW MC VAL: "<<binC <<std::endl;
165  return aux_par;
166  }
167 
168  float K_data_cv = -1.0;
169  float K_data_sys = -1.0;
170  float K_data_sta = -1.0;
171  float K_data = -1.0;
172  float K_aux = -1.0;
173 
174  if(tar.Tar_pdg == 321){
175  K_data_cv = vbin_datacv_pip[pip_bin] *vbin_datacv_kap_pip[binID];
176  K_data_sys = vbin_datasys_pip[pip_bin]*vbin_datasys_kap_pip[binID];
177  K_data_sta = vbin_datasta_pip[pip_bin]*vbin_datasta_kap_pip[binID];
178  K_data = K_data_sys + K_data_sta - K_data_cv;
179 
180  }
181  else if(tar.Tar_pdg ==-321){
182  K_data_cv = vbin_datacv_pim[pim_bin] *vbin_datacv_kam_pim[binID];
183  K_data_sys = vbin_datasys_pim[pim_bin]*vbin_datasys_kam_pim[binID];
184  K_data_sta = vbin_datasta_pim[pim_bin]*vbin_datasta_kam_pim[binID];
185  K_data = K_data_sys + K_data_sta - K_data_cv;
186  }
187  else if(tar.Tar_pdg ==310 || tar.Tar_pdg ==130){
188  //pip:
189  K_data_cv = vbin_datacv_pip[pip_bin] *vbin_datacv_kap_pip[binID];
190  K_data_sys = vbin_datasys_pip[pip_bin]*vbin_datasys_kap_pip[binID];
191  K_data_sta = vbin_datasta_pip[pip_bin]*vbin_datasta_kap_pip[binID];
192  K_aux = K_data_sys + K_data_sta - K_data_cv;
193  //pim:
194  K_data_cv = vbin_datacv_pim[pim_bin] *vbin_datacv_kam_pim[binID];
195  K_data_sys = vbin_datasys_pim[pim_bin]*vbin_datasys_kam_pim[binID];
196  K_data_sta = vbin_datasta_pim[pim_bin]*vbin_datasta_kam_pim[binID];
197  K_data = K_data_sys + K_data_sta - K_data_cv;
198  K_data *= 3.;
199  K_data += K_aux;
200  K_data *= 0.25;
201  }
202 
203  //check:
204  if(K_data_cv<low_value || K_data_sys<low_value || K_data_sta<low_value || K_data<low_value){
205  return 1.0;
206  }
207 
208  wgt = double(K_data)/binC;
209 
210  if(wgt<low_value){
211  // std::cout<<"TTMIPPK check wgt(<0) "<<iUniv<<" "<<tar.Pz<<" "<<tar.Pt<<" "<<tar.Tar_pdg<<std::endl;
212  return aux_par;
213  }
214  return wgt;
215 
216  }
NeutrinoFluxReweight::MIPPNumiMC::getInstance
static MIPPNumiMC * getInstance()
Definition: MIPPNumiMC.cpp:295
NeutrinoFluxReweight::MIPPNumiYieldsBins::getInstance
static MIPPNumiYieldsBins * getInstance()
Definition: MIPPNumiYieldsBins.cpp:153
std::vector< bool > NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::canReweight ( const InteractionChainData &  aa)
virtual

Look through the InteractionChainData input and identify those Interactions that can be reweighted as part of a chain. We return a vector indicating which elements will be assigned a weight by calculateWeight.

Implements NeutrinoFluxReweight::IInteractionChainReweighting.

Definition at line 86 of file MIPPNumiKaonYieldsReweighter.cpp.

86  {
87 
88  MIPPNumiYieldsBins* MIPPbins = MIPPNumiYieldsBins::getInstance();
89  std::vector<bool> this_nodes;
90  for(size_t ii=0;ii<(aa.interaction_chain).size();ii++){
91  this_nodes.push_back(false);
92  }
93 
94  //Look for MIPP Numi events for kaons
95  //if the code find a MIPP Numi event, it will look
96  //for how many interaction nodes covers
97  //if not, return all nodes false.
98  //bool is_there_mipp = false;
99  TargetData tar = aa.tar_info;
100 
101  //Cheking if the particle is a kaon plus or kaon minus or neutral kaon:
102  if(tar.Tar_pdg != 321 && tar.Tar_pdg != -321 && tar.Tar_pdg != 130 && tar.Tar_pdg != 310)return this_nodes;
103 
104  //kinematic coverage:
105  if(tar.Pz<20.0 || tar.Pz>80.0 || tar.Pt>2.0)return this_nodes;
106 
107  //data:
108  //Kaons:
109  int binID = MIPPbins->BinID(tar.Pz,tar.Pt,tar.Tar_pdg);
110  if(binID<0) return this_nodes;
111 
112  //Looking for a pion data:
113  int pip_bin = MIPPbins->BinID(tar.Pz,tar.Pt, 211);
114  int pim_bin = MIPPbins->BinID(tar.Pz,tar.Pt,-211);
115  if(tar.Tar_pdg == 321 && pip_bin<0)return this_nodes;
116  if(tar.Tar_pdg ==-321 && pim_bin<0)return this_nodes;
117  if(tar.Tar_pdg == 130 || tar.Tar_pdg == 310){
118  if(pip_bin<0 || pim_bin<0)return this_nodes;
119  }
120 
121  //Now that we know that we have a MIPP Numi event,
122  //we will see how many nodes are covered.
123  std::vector<InteractionData> this_interactions = aa.interaction_chain;
124 
125  //Now we have the index of the hadron that exit the target in the
126  //ancesty chain:
127  if(tar.Idx_ancestry>=0){
128  for(int ii=0;ii<tar.Idx_ancestry;ii++){
129  this_nodes[ii] = true;
130  }
131  }
132 
133  return this_nodes;
134  }
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
NeutrinoFluxReweight::MIPPNumiYieldsBins::getInstance
static MIPPNumiYieldsBins * getInstance()
Definition: MIPPNumiYieldsBins.cpp:153

Member Data Documentation

float NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::aux_par
private

Definition at line 35 of file MIPPNumiKaonYieldsReweighter.h.

const ParameterTable& NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::cvPars

Definition at line 27 of file MIPPNumiKaonYieldsReweighter.h.

int NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::iUniv
private

Definition at line 31 of file MIPPNumiKaonYieldsReweighter.h.

float NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::prt_no_inter
private

Definition at line 32 of file MIPPNumiKaonYieldsReweighter.h.

const ParameterTable& NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::univPars

Definition at line 28 of file MIPPNumiKaonYieldsReweighter.h.

std::vector<float> NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::vbin_datacv_kam_pim
private

Definition at line 34 of file MIPPNumiKaonYieldsReweighter.h.

std::vector<float> NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::vbin_datacv_kap_pip
private

Definition at line 34 of file MIPPNumiKaonYieldsReweighter.h.

std::vector<float> NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::vbin_datacv_pim
private

Definition at line 33 of file MIPPNumiKaonYieldsReweighter.h.

std::vector<float> NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::vbin_datacv_pip
private

Definition at line 33 of file MIPPNumiKaonYieldsReweighter.h.

std::vector<float> NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::vbin_datasta_kam_pim
private

Definition at line 34 of file MIPPNumiKaonYieldsReweighter.h.

std::vector<float> NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::vbin_datasta_kap_pip
private

Definition at line 34 of file MIPPNumiKaonYieldsReweighter.h.

std::vector<float> NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::vbin_datasta_pim
private

Definition at line 33 of file MIPPNumiKaonYieldsReweighter.h.

std::vector<float> NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::vbin_datasta_pip
private

Definition at line 33 of file MIPPNumiKaonYieldsReweighter.h.

std::vector<float> NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::vbin_datasys_kam_pim
private

Definition at line 34 of file MIPPNumiKaonYieldsReweighter.h.

std::vector<float> NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::vbin_datasys_kap_pip
private

Definition at line 34 of file MIPPNumiKaonYieldsReweighter.h.

std::vector<float> NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::vbin_datasys_pim
private

Definition at line 33 of file MIPPNumiKaonYieldsReweighter.h.

std::vector<float> NeutrinoFluxReweight::MIPPNumiKaonYieldsReweighter::vbin_datasys_pip
private

Definition at line 33 of file MIPPNumiKaonYieldsReweighter.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11