ThinTargetMesonIncidentReweighter.cpp
Go to the documentation of this file.
1 
2 #include "ThinTargetMesonIncidentReweighter.h"
3 #include <iostream>
4 #include "ThinTargetBins.h"
5 
6 namespace NeutrinoFluxReweight{
7 
8  ThinTargetMesonIncidentReweighter::ThinTargetMesonIncidentReweighter(int iuniv, const ParameterTable& cv_pars, const ParameterTable& univ_pars):iUniv(iuniv),cvPars(cv_pars),univPars(univ_pars){
9 
10  ThinTargetBins* Thinbins = ThinTargetBins::getInstance();
11 
12  vbin_pip_inc_pip.reserve(Thinbins->GetNbins_meson_incident());
13  vbin_pip_inc_pim.reserve(Thinbins->GetNbins_meson_incident());
14  vbin_pip_inc_kap.reserve(Thinbins->GetNbins_meson_incident());
15  vbin_pip_inc_kam.reserve(Thinbins->GetNbins_meson_incident());
16  vbin_pip_inc_k0.reserve(Thinbins->GetNbins_meson_incident());
17  vbin_pip_inc_p.reserve(Thinbins->GetNbins_meson_incident());
18  vbin_pip_inc_n.reserve(Thinbins->GetNbins_meson_incident());
19  vbin_pim_inc_pip.reserve(Thinbins->GetNbins_meson_incident());
20  vbin_pim_inc_pim.reserve(Thinbins->GetNbins_meson_incident());
21  vbin_pim_inc_kap.reserve(Thinbins->GetNbins_meson_incident());
22  vbin_pim_inc_kam.reserve(Thinbins->GetNbins_meson_incident());
23  vbin_pim_inc_k0.reserve(Thinbins->GetNbins_meson_incident());
24  vbin_pim_inc_p.reserve(Thinbins->GetNbins_meson_incident());
25  vbin_pim_inc_n.reserve(Thinbins->GetNbins_meson_incident());
26  vbin_kap_inc_pip.reserve(Thinbins->GetNbins_meson_incident());
27  vbin_kap_inc_pim.reserve(Thinbins->GetNbins_meson_incident());
28  vbin_kap_inc_kap.reserve(Thinbins->GetNbins_meson_incident());
29  vbin_kap_inc_kam.reserve(Thinbins->GetNbins_meson_incident());
30  vbin_kap_inc_k0.reserve(Thinbins->GetNbins_meson_incident());
31  vbin_kap_inc_p.reserve(Thinbins->GetNbins_meson_incident());
32  vbin_kap_inc_n.reserve(Thinbins->GetNbins_meson_incident());
33  vbin_kam_inc_pip.reserve(Thinbins->GetNbins_meson_incident());
34  vbin_kam_inc_pim.reserve(Thinbins->GetNbins_meson_incident());
35  vbin_kam_inc_kap.reserve(Thinbins->GetNbins_meson_incident());
36  vbin_kam_inc_kam.reserve(Thinbins->GetNbins_meson_incident());
37  vbin_kam_inc_k0.reserve(Thinbins->GetNbins_meson_incident());
38  vbin_kam_inc_p.reserve(Thinbins->GetNbins_meson_incident());
39  vbin_kam_inc_n.reserve(Thinbins->GetNbins_meson_incident());
40  vbin_k0_inc_pip.reserve(Thinbins->GetNbins_meson_incident());
41  vbin_k0_inc_pim.reserve(Thinbins->GetNbins_meson_incident());
42  vbin_k0_inc_kap.reserve(Thinbins->GetNbins_meson_incident());
43  vbin_k0_inc_kam.reserve(Thinbins->GetNbins_meson_incident());
44  vbin_k0_inc_k0.reserve(Thinbins->GetNbins_meson_incident());
45  vbin_k0_inc_p.reserve(Thinbins->GetNbins_meson_incident());
46  vbin_k0_inc_n.reserve(Thinbins->GetNbins_meson_incident());
47 
48  // const boost::interprocess::flat_map<std::string, double>& cv_table = cvPars.getMap();
49  // const boost::interprocess::flat_map<std::string, double>& univ_table = univPars.getMap();
50  char namepar[100];
51 
52  //meson left over:
53  sprintf(namepar,"ThinTarget_mesonleftover_incident_%d",0);
54  double dataval = univPars.getParameterValue(std::string(namepar));
55  bin_mesonleftover_inc = dataval;
56 
57  //5 incident particles, 7 produced particles:
58  const char* cinc[5] = {"pip","pim","kap","kam","k0"};
59  const char* cpro[7] = {"pip","pim","kap","kam","k0","p","n"};
60 
61  for(int ii=0;ii<5;ii++){
62  for(int jj=0;jj<7;jj++){
63  for(int kk=0;kk<4;kk++){
64 
65  sprintf(namepar,"ThinTarget_%s_incident_%s_%d",cinc[ii],cpro[jj],kk);
66  dataval = univPars.getParameterValue(std::string(namepar));
67  if(ii==0 && jj==0)vbin_pip_inc_pip.push_back(dataval);
68  if(ii==0 && jj==1)vbin_pip_inc_pim.push_back(dataval);
69  if(ii==0 && jj==2)vbin_pip_inc_kap.push_back(dataval);
70  if(ii==0 && jj==3)vbin_pip_inc_kam.push_back(dataval);
71  if(ii==0 && jj==4)vbin_pip_inc_k0.push_back(dataval);
72  if(ii==0 && jj==5)vbin_pip_inc_p.push_back(dataval);
73  if(ii==0 && jj==6)vbin_pip_inc_n.push_back(dataval);
74  if(ii==1 && jj==0)vbin_pim_inc_pip.push_back(dataval);
75  if(ii==1 && jj==1)vbin_pim_inc_pim.push_back(dataval);
76  if(ii==1 && jj==2)vbin_pim_inc_kap.push_back(dataval);
77  if(ii==1 && jj==3)vbin_pim_inc_kam.push_back(dataval);
78  if(ii==1 && jj==4)vbin_pim_inc_k0.push_back(dataval);
79  if(ii==1 && jj==5)vbin_pim_inc_p.push_back(dataval);
80  if(ii==1 && jj==6)vbin_pim_inc_n.push_back(dataval);
81  if(ii==2 && jj==0)vbin_kap_inc_pip.push_back(dataval);
82  if(ii==2 && jj==1)vbin_kap_inc_pim.push_back(dataval);
83  if(ii==2 && jj==2)vbin_kap_inc_kap.push_back(dataval);
84  if(ii==2 && jj==3)vbin_kap_inc_kam.push_back(dataval);
85  if(ii==2 && jj==4)vbin_kap_inc_k0.push_back(dataval);
86  if(ii==2 && jj==5)vbin_kap_inc_p.push_back(dataval);
87  if(ii==2 && jj==6)vbin_kap_inc_n.push_back(dataval);
88  if(ii==3 && jj==0)vbin_kam_inc_pip.push_back(dataval);
89  if(ii==3 && jj==1)vbin_kam_inc_pim.push_back(dataval);
90  if(ii==3 && jj==2)vbin_kam_inc_kap.push_back(dataval);
91  if(ii==3 && jj==3)vbin_kam_inc_kam.push_back(dataval);
92  if(ii==3 && jj==4)vbin_kam_inc_k0.push_back(dataval);
93  if(ii==3 && jj==5)vbin_kam_inc_p.push_back(dataval);
94  if(ii==3 && jj==6)vbin_kam_inc_n.push_back(dataval);
95  if(ii==4 && jj==0)vbin_k0_inc_pip.push_back(dataval);
96  if(ii==4 && jj==1)vbin_k0_inc_pim.push_back(dataval);
97  if(ii==4 && jj==2)vbin_k0_inc_kap.push_back(dataval);
98  if(ii==4 && jj==3)vbin_k0_inc_kam.push_back(dataval);
99  if(ii==4 && jj==4)vbin_k0_inc_k0.push_back(dataval);
100  if(ii==4 && jj==5)vbin_k0_inc_p.push_back(dataval);
101  if(ii==4 && jj==6)vbin_k0_inc_n.push_back(dataval);
102  }
103  }
104  }
105  }
106 
107  ThinTargetMesonIncidentReweighter::~ThinTargetMesonIncidentReweighter(){
108 
109  }
110  bool ThinTargetMesonIncidentReweighter::canReweight(const InteractionData& aa){
111  /*
112  if(aa.Inc_pdg != 211 && aa.Inc_pdg != -211 && aa.Inc_pdg != 321 && aa.Inc_pdg != -321 && aa.Inc_pdg != 130 && aa.Inc_pdg != 310)return false;
113  if(aa.Prod_pdg != 211 && aa.Prod_pdg != -211 && aa.Prod_pdg != 321 && aa.Prod_pdg != -321 && aa.Prod_pdg != 130 && aa.Prod_pdg != 310 && aa.Prod_pdg != 2212 && aa.Prod_pdg != 2112)return false;
114  */
115  if(aa.Proc.find("Inelastic")>100)return false;
116 
117  // ThinTargetBins* Thinbins = ThinTargetBins::getInstance();
118  //int bin = Thinbins->meson_inc_BinID(aa.xF,aa.Pt,aa.Prod_pdg);
119 
120  bool is_mesoninc = (aa.Inc_pdg >99 && aa.Inc_pdg < 1000) || (aa.Inc_pdg <-99 && aa.Inc_pdg > -1000);
121 
122  // if(bin>=0 || is_mesoninc)return true;
123  if(is_mesoninc)return true;
124  else return false;
125 
126  }
127 
128  double ThinTargetMesonIncidentReweighter::calculateWeight(const InteractionData& aa){
129 
130  double wgt = 1.0;
131 
132  bool right_inc = aa.Inc_pdg == 211 || aa.Inc_pdg == -211 || aa.Inc_pdg == 321 || aa.Inc_pdg == -321 || aa.Inc_pdg == 130 || aa.Inc_pdg == 310;
133  bool right_prod = aa.Prod_pdg == 211 || aa.Prod_pdg == -211 || aa.Prod_pdg == 321 || aa.Prod_pdg == -321 || aa.Prod_pdg == 130 || aa.Prod_pdg == 310 || aa.Prod_pdg == 2212 || aa.Prod_pdg == 2112;
134  ThinTargetBins* Thinbins = ThinTargetBins::getInstance();
135  int bin = Thinbins->meson_inc_BinID(aa.xF,aa.Pt,aa.Prod_pdg);
136  bool is_mesoninc = (aa.Inc_pdg >99 && aa.Inc_pdg < 1000) || (aa.Inc_pdg <-99 && aa.Inc_pdg > -1000);
137  if(bin>=0 && right_inc && right_prod){
138 
139  if(aa.Inc_pdg == 211){
140  if(aa.Prod_pdg == 211) wgt = vbin_pip_inc_pip[bin];
141  if(aa.Prod_pdg ==-211) wgt = vbin_pip_inc_pim[bin];
142  if(aa.Prod_pdg == 321) wgt = vbin_pip_inc_kap[bin];
143  if(aa.Prod_pdg ==-321) wgt = vbin_pip_inc_kam[bin];
144  if(aa.Prod_pdg ==130 || aa.Prod_pdg ==310) wgt = vbin_pip_inc_k0[bin];
145  if(aa.Prod_pdg ==2212) wgt = vbin_pip_inc_p[bin];
146  if(aa.Prod_pdg ==2112) wgt = vbin_pip_inc_n[bin];
147  }
148  else if(aa.Inc_pdg ==-211){
149  if(aa.Prod_pdg == 211) wgt = vbin_pim_inc_pip[bin];
150  if(aa.Prod_pdg ==-211) wgt = vbin_pim_inc_pim[bin];
151  if(aa.Prod_pdg == 321) wgt = vbin_pim_inc_kap[bin];
152  if(aa.Prod_pdg ==-321) wgt = vbin_pim_inc_kam[bin];
153  if(aa.Prod_pdg ==130 || aa.Prod_pdg ==310) wgt = vbin_pim_inc_k0[bin];
154  if(aa.Prod_pdg ==2212) wgt = vbin_pim_inc_p[bin];
155  if(aa.Prod_pdg ==2112) wgt = vbin_pim_inc_n[bin];
156  }
157  else if(aa.Inc_pdg == 321){
158  if(aa.Prod_pdg == 211) wgt = vbin_kap_inc_pip[bin];
159  if(aa.Prod_pdg ==-211) wgt = vbin_kap_inc_pim[bin];
160  if(aa.Prod_pdg == 321) wgt = vbin_kap_inc_kap[bin];
161  if(aa.Prod_pdg ==-321) wgt = vbin_kap_inc_kam[bin];
162  if(aa.Prod_pdg ==130 || aa.Prod_pdg ==310) wgt = vbin_kap_inc_k0[bin];
163  if(aa.Prod_pdg ==2212) wgt = vbin_kap_inc_p[bin];
164  if(aa.Prod_pdg ==2112) wgt = vbin_kap_inc_n[bin];
165  }
166  else if(aa.Inc_pdg ==-321){
167  if(aa.Prod_pdg == 211) wgt = vbin_kam_inc_pip[bin];
168  if(aa.Prod_pdg ==-211) wgt = vbin_kam_inc_pim[bin];
169  if(aa.Prod_pdg == 321) wgt = vbin_kam_inc_kap[bin];
170  if(aa.Prod_pdg ==-321) wgt = vbin_kam_inc_kam[bin];
171  if(aa.Prod_pdg ==130 || aa.Prod_pdg ==310) wgt = vbin_kam_inc_k0[bin];
172  if(aa.Prod_pdg ==2212) wgt = vbin_kam_inc_p[bin];
173  if(aa.Prod_pdg ==2112) wgt = vbin_kam_inc_n[bin];
174  }
175  else if(aa.Inc_pdg == 130 || aa.Inc_pdg == 310){
176  if(aa.Prod_pdg == 211) wgt = vbin_k0_inc_pip[bin];
177  if(aa.Prod_pdg ==-211) wgt = vbin_k0_inc_pim[bin];
178  if(aa.Prod_pdg == 321) wgt = vbin_k0_inc_kap[bin];
179  if(aa.Prod_pdg ==-321) wgt = vbin_k0_inc_kam[bin];
180  if(aa.Prod_pdg ==130 || aa.Prod_pdg ==310) wgt = vbin_k0_inc_k0[bin];
181  if(aa.Prod_pdg ==2212) wgt = vbin_k0_inc_p[bin];
182  if(aa.Prod_pdg ==2112) wgt = vbin_k0_inc_n[bin];
183  }
184  else{
185  std::cout<<"MESINC Something is wrong with pdg_inc: "<< aa.Inc_pdg <<" "<<aa.Prod_pdg<<std::endl;
186  return wgt;
187  }
188  }
189  else if(is_mesoninc){
190  wgt = bin_mesonleftover_inc;
191  }
192 
193  if(wgt<0){
194  //std::cout<<"TTMESONINC check wgt(<0) "<<iUniv<<" "<<aa.Inc_P<<" "<<aa.xF<<" "<<aa.Pt<<" "<<aa.Prod_pdg<<std::endl;
195  return 1.0;
196  }
197  if(wgt>10){
198  std::cout<<"BIG WGT IN TTMESONINC "<<iUniv<<" "<<wgt<<" "<<aa.Inc_P<<" "<<aa.xF<<" "<<aa.Pt<<" "<<aa.Prod_pdg<<std::endl;
199  return 1.0;
200  }
201  return wgt;
202 
203  }
204 
205 
206 }
NeutrinoFluxReweight::ThinTargetBins
A class to manage the bin definitions for MIPP Numi Yields.
Definition: ThinTargetBins.h:13
NeutrinoFluxReweight::ParameterTable
A list/table of parameter names and values.
Definition: ParameterTable.h:30
string
std::string string
Definition: nybbler.cc:12
NeutrinoFluxReweight::InteractionData::Prod_pdg
int Prod_pdg
pdg code of the produced particle
Definition: InteractionData.h:38
NeutrinoFluxReweight::InteractionData::Pt
double Pt
Transversal momentum (GeV/c) of the produced particle.
Definition: InteractionData.h:71
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
NeutrinoFluxReweight::InteractionData
The information about a hadronic interaction needed to calculate weights.
Definition: InteractionData.h:19
NeutrinoFluxReweight::ThinTargetMesonIncidentReweighter::calculateWeight
virtual double calculateWeight(const InteractionData &aa)
calculate a weight for this interaction given the central value parameters and the parameters for thi...
Definition: ThinTargetMesonIncidentReweighter.cpp:128
NeutrinoFluxReweight::InteractionData::xF
double xF
Feynmann-x of the produced particle: .
Definition: InteractionData.h:62
NeutrinoFluxReweight::ThinTargetBins::meson_inc_BinID
int meson_inc_BinID(double xf, double pt, int pdgcode)
Return Pion incident bin.
Definition: ThinTargetBins.cpp:337
NeutrinoFluxReweight::InteractionData::Proc
std::string Proc
Interaction process.
Definition: InteractionData.h:86
NeutrinoFluxReweight::InteractionData::Inc_P
double Inc_P
Momentum magnitude of the incident particle.
Definition: InteractionData.h:41
NeutrinoFluxReweight::ThinTargetMesonIncidentReweighter::ThinTargetMesonIncidentReweighter
ThinTargetMesonIncidentReweighter(int iuniv, const ParameterTable &cv_pars, const ParameterTable &univ_pars)
Definition: ThinTargetMesonIncidentReweighter.cpp:8
NeutrinoFluxReweight::InteractionData::Inc_pdg
int Inc_pdg
pdg code of the incident particle
Definition: InteractionData.h:35
NeutrinoFluxReweight::ThinTargetBins::getInstance
static ThinTargetBins * getInstance()
Definition: ThinTargetBins.cpp:406
bin
QTextStream & bin(QTextStream &s)
Definition: qtextstream.cpp:2006
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
NeutrinoFluxReweight::ThinTargetMesonIncidentReweighter::canReweight
virtual bool canReweight(const InteractionData &aa)
can the particular instance of this class reweight this interaction?
Definition: ThinTargetMesonIncidentReweighter.cpp:110
Generated by   doxygen 1.8.11