All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Attributes | List of all members
sbn::evwgh::FluxWeightCalc Class Reference

#include <FluxCalcPrep.h>

Inheritance diagram for sbn::evwgh::FluxWeightCalc:
sbn::evwgh::WeightCalc

Public Member Functions

 FluxWeightCalc ()
 
void Configure (fhicl::ParameterSet const &p, CLHEP::HepRandomEngine &engine) override
 
std::vector< float > GetWeight (art::Event &e, size_t inu) override
 
double UnisimWeightCalc (double enu, int ptype, int ntype, double randomN, bool noNeg)
 
std::pair< bool, double > PHNWeightCalc (simb::MCFlux flux, float rand)
 
std::pair< bool, double > PHFSWeightCalc (simb::MCFlux flux, std::vector< float > rand)
 
std::pair< bool, double > PHSWWeightCalc (simb::MCFlux flux, std::vector< float > rand)
 
std::pair< bool, double > PHSWCSVWeightCalc (simb::MCFlux flux, std::vector< float > rand)
 
std::vector< double > ConvertToVector (TArrayD const *array)
 
- Public Member Functions inherited from sbn::evwgh::WeightCalc
void SetName (std::string name)
 
void SetType (std::string type)
 
std::string GetName ()
 
std::string GetType ()
 
std::string GetFullName ()
 

Private Attributes

std::string fGeneratorModuleLabel
 
std::string CalcType
 
double fScalePos {}
 
double fScaleNeg = 1
 
double fCV [4][4][200]
 
double fRWpos [4][4][200]
 
double fRWneg [4][4][200]
 
bool PosOnly {false}
 
std::vector< int > fprimaryHad = {0}
 
std::vector< double > FitVal {}
 
TMatrixD * FitCov {nullptr}
 
TMatrixD * HARPXSec {nullptr}
 
std::vector< double > HARPmomentumBounds {}
 
std::vector< double > HARPthetaBounds {}
 
std::vector< double > SWParam {}
 
bool fIsDecomposed {false}
 
int wc = 0
 
int wcn = 0
 
int wc0 = 0
 
int wc1 = 0
 
int wc30 = 0
 

Additional Inherited Members

- Public Attributes inherited from sbn::evwgh::WeightCalc
EventWeightParameterSet fParameterSet
 

Detailed Description

Definition at line 22 of file FluxCalcPrep.h.

Constructor & Destructor Documentation

sbn::evwgh::FluxWeightCalc::FluxWeightCalc ( )
inline

Definition at line 24 of file FluxCalcPrep.h.

24 : WeightCalc() {}

Member Function Documentation

void sbn::evwgh::FluxWeightCalc::Configure ( fhicl::ParameterSet const &  p,
CLHEP::HepRandomEngine &  engine 
)
overridevirtual

Implements sbn::evwgh::WeightCalc.

Definition at line 11 of file FluxCalcPrep.cxx.

12  {
13  std::cout<<"SBNEventWeight Flux: configure "<< GetName() << std::endl;
14 
15  fGeneratorModuleLabel = p.get<std::string>("generator_module_label");//use this label to get MC*Handle
16  const fhicl::ParameterSet& pset = p.get<fhicl::ParameterSet>(GetName());
17 
18  //0. << Reweighting Environment >>
19 
20  //skip "random_seed"
21  auto const& pars = pset.get<std::vector<std::string> >("parameter_list");
22  std::vector< float > parsigmas(pars.size(), 1.0);
23  if (pars.size() != parsigmas.size()) {
24  throw cet::exception(__PRETTY_FUNCTION__) << GetName() << ": "
25  << "parameter_list and parameter_sigma length mismatch."
26  << std::endl;
27  }
28 
29  if(!pset.get_if_present("parameter_sigma", parsigmas)){
30  std::cout<<"SBNEventWeight Flux: `parameter_sigma` was not set; now it is set as 1"<<std::endl;
31  }
32 
33  std::string fMode = pset.get<std::string>("mode");//3 types: multisim/pmNsigma/fixed
34  int number_of_multisims = pset.get<int>("number_of_multisims", 1);
35 
36  fParameterSet.Configure(GetFullName(), fMode, number_of_multisims);
37  for (size_t i=0; i<pars.size(); i++) {
38  //Note: no sigma is used in this script.
39  fParameterSet.AddParameter(pars[i], parsigmas[i]);//parsigmas[i]);
40  }
41 
42 
43  //1. << Calculator Settings >>
44  CalcType = pset.get< std::string >("calc_type");//Unisim,PrimaryHadronSWCentralSplineVariation,PrimaryHadronFeynmanScaling,PrimaryHadronSanfordWang,PrimaryHadronNormalization
45  std::cout<<"SBNEventWeight Flux: Calculator type: "<<CalcType<<std::endl;
46 
47  fScalePos = pset.get<double>("scale_factor_pos");
48  if(!pset.get_if_present("scale_factor_neg",fScaleNeg)){
49  std::cout<<"SBNEventWeight Flux: auto-assignment: scale_factor_neg = 1."<<std::endl;
50  }
51 
52  cet::search_path sp("FW_SEARCH_PATH");
53  bool validC = true;
54 
55  //----------------------------
56  //-- Non hadrons production --
57  //----------------------------
58 
59 
60  fParameterSet.Sample(engine);//random_seed is loaded at sbncode/Base/WeightManager.h
61 
62  if( CalcType == "Unisim"){//Unisim Calculator
63 
64  std::string dataInput1 = pset.get< std::string >("CentralValue_hist_file");
65  std::string cvfile_path = sp.find_file(dataInput1);
66  TFile fcv(Form("%s",cvfile_path.c_str()));
67 
68  std::string dataInput2pos = pset.get< std::string >("PositiveSystematicVariation_hist_file");
69  std::string rwfilepos = sp.find_file(dataInput2pos);
70  TFile frwpos(Form("%s", rwfilepos.c_str()));
71 
72  std::string dataInput2neg = pset.get< std::string >("NegativeSystematicVariation_hist_file");
73  std::string rwfileneg = sp.find_file(dataInput2neg);
74  TFile frwneg(Form("%s", rwfileneg.c_str()));
75 
76  if(dataInput2pos == dataInput2neg) PosOnly = true;//true - for skin depth
77  //Those May07*.root use the following convention to name histograms
78  int cptype[4] = {1,2,3,4}; //mu, pi, k0, k
79  int cntype[4] = {1,2,3,4}; //nue, anue, numu, anumu
80 
81  for (int iptyp=0;iptyp<4;iptyp++) {
82  for (int intyp=0;intyp<4;intyp++) {
83  for (int ibin=0;ibin<200;ibin++) { //Grab events from ibin+1
84  fCV[iptyp][intyp][ibin]=(dynamic_cast<TH1F*> (fcv.Get(Form("h5%d%d",cptype[iptyp],cntype[intyp]))))->GetBinContent(ibin+1);
85  fRWpos[iptyp][intyp][ibin]=(dynamic_cast<TH1F*> (frwpos.Get(Form("h5%d%d",cptype[iptyp],cntype[intyp]))))->GetBinContent(ibin+1);
86  fRWneg[iptyp][intyp][ibin]=(dynamic_cast<TH1F*> (frwneg.Get(Form("h5%d%d",cptype[iptyp],cntype[intyp]))))->GetBinContent(ibin+1);
87  }// energy bin
88  }// type of neutrinos
89  }//type of hadron parent
90  fcv.Close();
91  frwpos.Close();
92  frwneg.Close();
93 
94  //-------------
95  //-- Hadrons --
96  //-------------
97  } else if( CalcType.compare(0, 13,"PrimaryHadron") == 0){//Hadron Calculators
98 
99  fprimaryHad = pset.get< std::vector<int>>("PrimaryHadronGeantCode");
100 
101  if( CalcType == "PrimaryHadronNormalization" ){//k-
102 
103  fParameterSet.Sample(engine);//Yes.. load it again;
104 
105  } else{//Other Hadron Calculators
106 
107  std::string dataInput = pset.get< std::string >("ExternalData");
108  std::string ExternalDataInput = sp.find_file(dataInput);
109  TFile* file = new TFile(Form("%s",ExternalDataInput.c_str()));
110 
111  std::vector< std::string > pname; // these are keys to histograms
112  if( CalcType == "PrimaryHadronFeynmanScaling" ){//k+
113 
114  pname.push_back("FS/KPlus/FSKPlusFitVal");
115  pname.push_back("FS/KPlus/FSKPlusFitCov");
116 
117  TArrayD* FSKPlusFitValArray = (TArrayD*) file->Get(pname[0].c_str());
118  FitVal = FluxWeightCalc::ConvertToVector(FSKPlusFitValArray);
119  FitCov = (TMatrixD*) file->Get(pname[1].c_str());
120  *(FitCov) *= fScalePos*fScalePos;
121 
122  }else if( CalcType == "PrimaryHadronSanfordWang" ){//k0
123  pname.push_back("SW/K0s/SWK0sFitVal");
124  pname.push_back("SW/K0s/SWK0sFitCov");
125  TArrayD* SWK0FitValArray = (TArrayD*) file->Get(pname[0].c_str());
126 
127  FitVal = FluxWeightCalc::ConvertToVector(SWK0FitValArray);//TArrayD--> vector
128  FitCov = (TMatrixD*) file->Get(pname[1].c_str());
129  *(FitCov) *= fScalePos*fScalePos;
130 
131  }else if( CalcType == "PrimaryHadronSWCentralSplineVariation" ){//pi+-
132 
133  std::string fitInput = pset.get< std::string >("ExternalFit");
134  std::string HadronName;
135  std::string HadronAbriviation;
136  if(fprimaryHad[0] == 211){
137  HadronName = "PiPlus";
138  HadronAbriviation = "PP";
139  } else if(fprimaryHad[0] == -211){
140  HadronName = "PiMinus";
141  HadronAbriviation = "PM";
142  } else{
143  throw art::Exception(art::errors::StdException)
144  << "sanford-wang is only configured for charged pions ";
145  }
146  pname.push_back( Form("HARPData/%s/%sCrossSection",HadronName.c_str(),HadronAbriviation.c_str()) ); // Cross Section
147  pname.push_back( Form("HARPData/%s/%scovarianceMatrix",HadronName.c_str(),HadronAbriviation.c_str()) ); // Covariance Matrix
148  pname.push_back( Form("HARPData/%s/%smomentumBoundsArray",HadronName.c_str(),HadronAbriviation.c_str()) ); // Momentum Bounds
149  pname.push_back( Form("HARPData/%s/%sthetaBoundsArray",HadronName.c_str(),HadronAbriviation.c_str()) ); // Theta Bounds
150 
151  HARPXSec = (TMatrixD*) file->Get(pname[0].c_str());
152  TMatrixD* HARPCov = (TMatrixD*) file->Get(pname[1].c_str());
153 
154  TDecompChol dc = TDecompChol(*(HARPCov));//perform Choleskey Decomposition
155  if(!dc.Decompose()){
156  throw art::Exception(art::errors::StdException)
157  << "Cannot decompose covariance matrix to begin smearing.";
158  }
159 
160  //Get upper triangular matrix. This maintains the relations in the
161  // covariance matrix, but simplifies the structure.
162  fIsDecomposed = true;
163  FitCov = new TMatrixD(dc.GetU());
164 
165 
166  TArrayD* HARPmomentumBoundsArray = (TArrayD*) file->Get(pname[2].c_str());
167  HARPmomentumBounds = FluxWeightCalc::ConvertToVector(HARPmomentumBoundsArray);
168 
169  TArrayD* HARPthetaBoundsArray = (TArrayD*) file->Get(pname[3].c_str());
170  HARPthetaBounds = FluxWeightCalc::ConvertToVector(HARPthetaBoundsArray);
171 
172  /////////////////
173  //
174  // Extract the Sanford-Wang Fit Parmeters
175  //
176  ////////////////
177 
178  std::string ExternalFitInput = sp.find_file(fitInput);
179  TFile* Fitfile = new TFile(Form("%s",ExternalFitInput.c_str()));
180 
181  std::string fitname; // these are what we will extract from the file
182  fitname = Form("SW/%s/SW%sFitVal",HadronName.c_str(),HadronName.c_str()); // Sanford-Wang Fit Parameters
183 
184  TArrayD* SWParamArray = (TArrayD*) Fitfile->Get(fitname.c_str());
185  SWParam = FluxWeightCalc::ConvertToVector(SWParamArray);
186 
187 
188  }else validC = false;//slightly incorrect calculator name in *fcl
189 
190  if(validC){//load random numbers based on FitCov
191  //{2*multisim vector}
192  //{{Ncols() elements}} <-- feed these amount everytime;
193  for( int index = 1; index < 2*(FitCov->GetNcols()); index ++){
194  fParameterSet.Sample(engine);//load it 2*<number_of_multisims> times
195  }
196 
197  }else {
198  throw cet::exception(__PRETTY_FUNCTION__) << GetName() << ": "
199  <<" calculator "+CalcType + "is invalid"
200  <<std::endl;
201  }
202 
203 
204  }//end of special Hadron calculator configurations
205 
206  } else validC = false; //the calculator name is way too off.
207 // std::cout<<"SBNEventWeight : finish configuration."<<std::endl;
208  }//End of Configure() function
sbn::evwgh::EventWeightParameterSet::AddParameter
void AddParameter(std::string name, float width, float mean=0, size_t covIndex=0)
Definition: EventWeightParameterSet.cxx:51
sbn::evwgh::FluxWeightCalc::FitVal
std::vector< double > FitVal
Definition: FluxCalcPrep.h:71
sbn::evwgh::FluxWeightCalc::HARPmomentumBounds
std::vector< double > HARPmomentumBounds
Definition: FluxCalcPrep.h:76
sbn::evwgh::FluxWeightCalc::fGeneratorModuleLabel
std::string fGeneratorModuleLabel
Definition: FluxCalcPrep.h:55
file
* file
Definition: file_to_url.sh:69
p
pdgs p
Definition: selectors.fcl:22
sbn::evwgh::FluxWeightCalc::fRWpos
double fRWpos[4][4][200]
Definition: FluxCalcPrep.h:64
sbn::evwgh::EventWeightParameterSet::Configure
void Configure(std::string name, ReweightType rwtype, size_t nuni=1)
Definition: EventWeightParameterSet.cxx:12
evgen::fMode
createEngine fMode
Definition: GaisserParam_module.cc:168
sbn::evwgh::EventWeightParameterSet::Sample
void Sample(CLHEP::HepRandomEngine &engine)
Definition: EventWeightParameterSet.cxx:63
sbn::evwgh::FluxWeightCalc::SWParam
std::vector< double > SWParam
Definition: FluxCalcPrep.h:78
sbn::evwgh::FluxWeightCalc::ConvertToVector
std::vector< double > ConvertToVector(TArrayD const *array)
Definition: FluxCalcPrep.cxx:354
sbn::evwgh::FluxWeightCalc::HARPthetaBounds
std::vector< double > HARPthetaBounds
Definition: FluxCalcPrep.h:77
sbn::evwgh::FluxWeightCalc::fRWneg
double fRWneg[4][4][200]
Definition: FluxCalcPrep.h:65
sbn::evwgh::FluxWeightCalc::fprimaryHad
std::vector< int > fprimaryHad
Definition: FluxCalcPrep.h:69
sbn::evwgh::FluxWeightCalc::fCV
double fCV[4][4][200]
Definition: FluxCalcPrep.h:63
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
std::vector< double > sbn::evwgh::FluxWeightCalc::ConvertToVector ( TArrayD const *  array)

Definition at line 354 of file FluxCalcPrep.cxx.

354  {
355  std::vector<double> v(array->GetSize());
356  std::copy(array->GetArray(), array->GetArray() + array->GetSize(),
357  v.begin());
358  return v;
359  } // ConvertToVector()
array
then echo ***************************************echo array
Definition: find_fhicl.sh:28
copy
T copy(T const &v)
Definition: RangeForWrapper_test.cc:98
std::vector< float > sbn::evwgh::FluxWeightCalc::GetWeight ( art::Event &  e,
size_t  inu 
)
overridevirtual

Implements sbn::evwgh::WeightCalc.

Definition at line 211 of file FluxCalcPrep.cxx.

211  {
212  bool count_weights = false;
213 // std::cout<<"SBNEventWeight : getweight for the "<<inu<<" th particles of an event"<< std::endl;
214  //MCFlux & MCTruth
215  art::Handle< std::vector<simb::MCFlux> > mcFluxHandle;
216  e.getByLabel(fGeneratorModuleLabel,mcFluxHandle);
217  std::vector<simb::MCFlux> const& fluxlist = *mcFluxHandle;
218  //or do the above 3 lines in one line
219 // auto const& mclist = *e.getValidHandle<std::vector<simb::MCTruth>>(fGeneratorModuleLabel);
220 
221  // If no neutrinos in this event, gives 0 weight;
222  int NUni = fParameterSet.fNuniverses;
223  std::vector<float> weights;//( mclist.size(), 0);
224  if(fluxlist.size() == 0){
225  std::cout<<"SBNEventWeight Flux: EMPTY WEIGHTS"<<std::endl;
226  return weights;
227  }
228 
229  //Iterate through each neutrino in the event
230  std::cout<<"SBNEventWeight Flux: Get calculator "<<CalcType<<std::endl;
231  if( CalcType == "Unisim"){//Unisim Calculator
232  weights.resize(NUni);
233  //Unisim specific
234  //containers for the parent and neutrino type information
235  int ptype = std::numeric_limits<int>::max();
236  int ntype = std::numeric_limits<int>::max();
237 
238  // Discover the neutrino parent type
239  // This contains the neutrino's parentage information
240  if ( fluxlist[inu].fptype==13 || fluxlist[inu].fptype==-13 ) ptype = 0;//mu
241  else if ( fluxlist[inu].fptype==211 || fluxlist[inu].fptype==-211 ) ptype = 1;//pi
242  else if ( fluxlist[inu].fptype==130 ) ptype = 2;//K0
243  else if ( fluxlist[inu].fptype==321 || fluxlist[inu].fptype==-321 ) ptype = 3;//K
244  else {
245  throw cet::exception(__FUNCTION__) << GetName()<<"::Unknown ptype "<<fluxlist[0].fptype<< std::endl;
246  }
247 
248  // Discover the neutrino type
249  // This contains the neutrino's flavor information
250  if ( fluxlist[inu].fntype== 12 ) ntype=0;//nue
251  else if ( fluxlist[inu].fntype==-12 ) ntype=1;//nuebar
252  else if ( fluxlist[inu].fntype== 14 ) ntype=2;//numu
253  else if ( fluxlist[inu].fntype==-14 ) ntype=3;//numubar
254  else {
255  throw cet::exception(__FUNCTION__) << GetName()<<"::Unknown ntype "<<fluxlist[inu].fntype<< std::endl;
256  }
257 
258  // Collect neutrino energy; mclist is replaced with fluxlist.
259 // double enu= mclist[inu].GetNeutrino().Nu().E();
260  double enu= fluxlist[inu].fnenergyn;
261 
262  //Let's make a weights based on the calculator you have requested
263 
264  if(fParameterSet.fRWType == EventWeightParameterSet::kMultisim){
265 
266  for (size_t i=0;i<weights.size();i++) {
267  double randomN = (fParameterSet.fParameterMap.begin())->second[i];
268  weights[i]=UnisimWeightCalc(enu, ptype, ntype, randomN , PosOnly);//AddParameter result does not work here;
269 
270  if(count_weights){
271  if(weights[i]<0){
272  wcn++;
273  }else if((weights[i]-0)<1e-30){
274  wc0++;
275  }else if(fabs(weights[i]-30) < 1e-30){
276  wc30++;
277  } else if(fabs(weights[i]-1)<1e-30){
278  wc1++;
279  } else {
280  wc++;
281  }
282  }
283  }//Iterate through the number of universes
284  }
285  } else{//then this must be PrimaryHadron
286 
287 
288  // First let's check that the parent of the neutrino we are looking for is
289  // the particle we intended it to be, if not set all weights to 1
290  if (std::find(fprimaryHad.begin(), fprimaryHad.end(),(fluxlist[inu].ftptype)) == fprimaryHad.end() ){//if it does not contain any particles we need get 1
291  weights.resize( NUni);
292  std::fill(weights.begin(), weights.end(), 1);
293  return weights;//done, all 1
294  }// Hadronic parent check
295 
296  if(fParameterSet.fRWType == EventWeightParameterSet::kMultisim){
297 
298  for (unsigned int i = 0; int(weights.size()) < NUni; i++) {//if all weights are 1, no need to calculate weights;
299  std::pair<bool, double> test_weight;
300 
301  std::vector< float > Vrandom = (fParameterSet.fParameterMap.begin())->second;//vector of random #
302  std::vector< float > subVrandom;//sub-vector of random numbers;
303  if( CalcType == "PrimaryHadronNormalization"){//Normalization
304  test_weight = PHNWeightCalc(fluxlist[inu], Vrandom[i]);
305 
306  } else {
307  subVrandom = {Vrandom.begin()+i*FitCov->GetNcols(), Vrandom.begin()+(i+1)*FitCov->GetNcols()};
308  if( CalcType == "PrimaryHadronFeynmanScaling"){//FeynmanScaling
309  test_weight = PHFSWeightCalc(fluxlist[inu], subVrandom);
310 
311  } else if( CalcType == "PrimaryHadronSanfordWang"){//SanfordWang
312  test_weight = PHSWWeightCalc(fluxlist[inu], subVrandom);
313 
314  } else if( CalcType == "PrimaryHadronSWCentralSplineVariation"){//SWCentaralSplineVariation
315  test_weight = PHSWCSVWeightCalc(fluxlist[inu], subVrandom);
316 
317  } else throw cet::exception(__PRETTY_FUNCTION__) << GetName() << ": this shouldnt happen.."<<std::endl;
318  }
319 
320  if(test_weight.first){
321  weights.push_back(test_weight.second);
322  double tmp_weight = test_weight.second;
323 // std::cout<<i<<" ("<<tmp_weight<<") ";
324  if(count_weights){
325  if(tmp_weight<0){
326  wcn++;
327  }else if((tmp_weight-0)<1e-30){
328  wc0++;
329  }else if(fabs(tmp_weight-30) < 1e-30){
330  wc30++;
331  } else if(fabs(tmp_weight-1)<1e-30){
332  wc1++;
333  } else {
334  wc++;
335  }
336  }
337 
338  };
339 
340  }//Iterate through the number of universes
341  }//Yes, Multisim
342  }
343 
344  if(count_weights){
345  std::cout<<"SBNEventWeight Flux: Weights counter: (normal, <0, ==0, 1, 30)= ";
346  std::cout<<wc<<", "<<wcn<<", "<<wc0<<", "<<wc1<<", "<<wc30<<std::endl;
347  }
348 // std::cout<<"Next partile/event\n"<<std::endl;
349 
350  return weights;
351  }//GetWeight()
second
second
Definition: detectorproperties_icarus.fcl:8
sbn::evwgh::FluxWeightCalc::fGeneratorModuleLabel
std::string fGeneratorModuleLabel
Definition: FluxCalcPrep.h:55
sbn::evwgh::FluxWeightCalc::PHSWCSVWeightCalc
std::pair< bool, double > PHSWCSVWeightCalc(simb::MCFlux flux, std::vector< float > rand)
Definition: PrimaryHadronSWCentralSplineVariationWeightCalc.cxx:10
sbn::evwgh::FluxWeightCalc::UnisimWeightCalc
double UnisimWeightCalc(double enu, int ptype, int ntype, double randomN, bool noNeg)
Definition: UnisimWeightCalc.cxx:9
sbn::evwgh::FluxWeightCalc::PHSWWeightCalc
std::pair< bool, double > PHSWWeightCalc(simb::MCFlux flux, std::vector< float > rand)
Definition: PrimaryHadronSanfordWangWeightCalc.cxx:9
sbn::evwgh::EventWeightParameterSet::fRWType
ReweightType fRWType
Type of throws (the same for all parameters in a set)
Definition: EventWeightParameterSet.h:146
sbn::evwgh::EventWeightParameterSet::fParameterMap
std::map< EventWeightParameter, std::vector< float > > fParameterMap
Mapping of definitions to the set of values.
Definition: EventWeightParameterSet.h:143
trkf::fill
void fill(const art::PtrVector< recob::Hit > &hits, int only_plane)
Definition: KHitContainer.cxx:18
sbn::evwgh::EventWeightParameterSet::fNuniverses
size_t fNuniverses
Number of universes (i.e. random throws)
Definition: EventWeightParameterSet.h:147
ntype
then ntype
Definition: condor_hadd_sam.sh:202
e
do i e
Definition: updateVersionList.sh:68
sbn::evwgh::FluxWeightCalc::PHNWeightCalc
std::pair< bool, double > PHNWeightCalc(simb::MCFlux flux, float rand)
Definition: PrimaryHadronNormalizationWeightCalc.cxx:9
sbn::evwgh::FluxWeightCalc::PHFSWeightCalc
std::pair< bool, double > PHFSWeightCalc(simb::MCFlux flux, std::vector< float > rand)
Definition: PrimaryHadronFeynmanScalingWeightCalc.cxx:9
sbn::evwgh::FluxWeightCalc::fprimaryHad
std::vector< int > fprimaryHad
Definition: FluxCalcPrep.h:69
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
std::pair< bool, double > sbn::evwgh::FluxWeightCalc::PHFSWeightCalc ( simb::MCFlux  flux,
std::vector< float >  rand 
)

Definition at line 9 of file PrimaryHadronFeynmanScalingWeightCalc.cxx.

9  {
10 
11  //
12  // Largely built off the MiniBooNE code
13  // but this is intended to expand beyond it
14  //
15  // Edits from MiniBooNE code:
16  //
17  // JZ (6/2017) : Remove max weight, set to max_limit of a double
18  // JZ (6/2017) : Fixed typo in the E_cm calculation per Mike S. suggestions
19  //
20 
21  // We need to guard against unphysical parameters
22  bool parameters_pass = false;
23 
24  // Get Neutrino Parent Kinimatics
25  double HadronMass = 0.4937;
26 
27  // if(fabs(fprimaryHad) == 321) HadronMass = 0.4937; //Charged Kaon
28  // else{
29  // throw art::Exception(art::errors::StdException)
30  // << "Feynman Scaling is only configured for Charged Kaons ";
31  // }
32 
33  TLorentzVector HadronVec;
34  double HadronPx = flux.ftpx;
35  double HadronPy = flux.ftpy;
36  double HadronPz = flux.ftpz;
37  double HadronE = sqrt(HadronPx*HadronPx +
38  HadronPy*HadronPy +
39  HadronPz*HadronPz +
40  HadronMass*HadronMass);
41  HadronVec.SetPxPyPzE(HadronPx,HadronPy,HadronPz,HadronE);
42 
43  // Get Initial Proton Kinitmatics
44  // CURRENTLY GSimple flux files drop information about
45  // the initial state proton, but that this
46  TLorentzVector ProtonVec;
47  double ProtonMass = 0.9382720;
48  double ProtonPx = 0;
49  double ProtonPy = 0;
50  double ProtonPz = 8.89; //GeV
51  double ProtonE = sqrt(ProtonPx*ProtonPx +
52  ProtonPy*ProtonPy +
53  ProtonPz*ProtonPz +
54  ProtonMass*ProtonMass);
55  ProtonVec.SetPxPyPzE(ProtonPx,ProtonPy,ProtonPz,ProtonE);
56 
57  //
58  // From Mike S. :
59  // The fitting was done without this error so I think that you need to use:
60  // ecm = sqrt(2.*mp**2+2.*mp*eb)
61  //
62  double E_cm = sqrt(2*ProtonMass*ProtonMass + 2*ProtonMass*ProtonVec.E());
63 
64  double gamma = (ProtonE + ProtonMass) / E_cm;
65  double gammaBeta = ProtonVec.P()/E_cm;
66 
67  double MaxThreshold = 2.053;
68  // if(fabs(fprimaryHad) == 321){
69  // MaxThreshold = 2.053; // Mike S. Comment: K+: lambda + proton
70  // }
71  // else{
72  // throw art::Exception(art::errors::StdException)
73  // << "Feynman Scaling is only configured for Charged Kaons ";
74  // }
75 
76  double E_cm_max = (E_cm*E_cm
77  +HadronMass*HadronMass
78  -MaxThreshold*MaxThreshold)/(2*E_cm);
79  double p_cm_max = sqrt(E_cm_max*E_cm_max-HadronMass*HadronMass);
80 
81  double HadronPT = HadronVec.Pt();
82  double HadronPparallel = HadronVec.P()*cos(HadronVec.Theta());
83  double ParallelCMp = gamma*HadronPparallel - gammaBeta*HadronE;
84 
85  double xF;
86 
87  xF = ParallelCMp / p_cm_max;
88 
89  // Define the scale central value cross section
90 
91  // Pull out the parameters for the Feynman Scaling
92  double c1 = FitVal.at(0);
93  double c2 = FitVal.at(1);
94  double c3 = FitVal.at(2);
95  double c4 = FitVal.at(3);
96  double c5 = FitVal.at(4);
97  double c6 = FitVal.at(5);
98  double c7 = FitVal.at(6);
99 
100  double CV = c1*(HadronVec.P()*HadronVec.P()/HadronE)*exp(-1.*c3*pow(fabs(xF),c4)
101  - c7*pow(fabs(HadronPT*xF),c6)
102  - c2*HadronPT
103  - c5*HadronPT*HadronPT);
104 
105  if(CV < 0) CV = 0;
106  if(fabs(xF) > 1) CV = 0;
107  // Define the variations around that cross section
108  // To do this we will call the a function from WeightCalc
109  // that will generate a set of smeared parameters based
110  // on the covariance matrix that we imported
111  std::vector< double > FSKPlusFitSmeared = MultiGaussianSmearing(FitVal, FitCov, rand);//Defined in sbncode/SBNEventWeight/Base/SmearingUtils.h
112 
113  // Pull out the parameters for the Feynman Scaling
114  double smeared_c1 = FSKPlusFitSmeared.at(0);
115  double smeared_c2 = FSKPlusFitSmeared.at(1);
116  double smeared_c3 = FSKPlusFitSmeared.at(2);
117  double smeared_c4 = FSKPlusFitSmeared.at(3);
118  double smeared_c5 = FSKPlusFitSmeared.at(4);
119  double smeared_c6 = FSKPlusFitSmeared.at(5);
120  double smeared_c7 = FSKPlusFitSmeared.at(6);
121 
122  if(smeared_c1 > 0 &&
123  smeared_c2 > 0 &&
124  smeared_c3 > 0 &&
125  smeared_c5 > 0 &&
126  smeared_c7 > 0){
127  parameters_pass = true;
128  }
129 
130  double RW = smeared_c1*(HadronVec.P()*HadronVec.P()/HadronE)*exp(-1.*smeared_c3*pow(fabs(xF),smeared_c4)
131  - smeared_c7*pow(fabs(HadronPT*xF),smeared_c6)
132  - smeared_c2*HadronPT
133  - smeared_c5*HadronPT*HadronPT);
134 
135  if(RW < 0) RW = 0;
136  if(fabs(xF) > 1) RW = 0;
137 
138  double weight = 1;
139 
140  if(RW < 0 || CV < 0){
141  weight = 1;
142  }
143  else if(CV < 1.e-12){
144  weight = 1;
145  }
146  else{
147  weight *= RW/CV;
148  }
149 
150  if(weight < 0) weight = 0;
151  if(weight > 30) weight = 30;
152  if(!(std::isfinite(weight))){
153  std::cout << "FS : Failed to get a finite weight" << std::endl;
154  weight = 30;}
155 
156  std::pair<bool, double> output(parameters_pass, weight);
157 
158 
159 
160  return output;
161 
162  }
sbn::evwgh::FluxWeightCalc::FitVal
std::vector< double > FitVal
Definition: FluxCalcPrep.h:71
geo::vect::isfinite
bool isfinite(Vector const &v)
Returns whether all components of the vector are finite.
Definition: geo_vectors_utils.h:2135
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
output
BEGIN_PROLOG sequence::SlidingWindowTriggerPatternsOppositeWindows END_PROLOG simSlidingORM6O6 effSlidingORW output
Definition: triggersim_eastmodule_icarus.fcl:134
sbn::evwgh::MultiGaussianSmearing
std::vector< std::vector< double > > MultiGaussianSmearing(std::vector< double > const &centralValue, std::vector< std::vector< double > > const &inputCovarianceMatrix, int n_multisims, CLHEP::RandGaussQ &GaussRandom)
Definition: SmearingUtils.cxx:10
e
do i e
Definition: updateVersionList.sh:68
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
std::pair< bool, double > sbn::evwgh::FluxWeightCalc::PHNWeightCalc ( simb::MCFlux  flux,
float  rand 
)

Definition at line 9 of file PrimaryHadronNormalizationWeightCalc.cxx.

9  {
10  // We need to guard against unphysical parameters
11  bool parameters_pass = false;
12 
13  double weight = rand+1;
14 
15  if(weight > 0) parameters_pass = true;
16  else{parameters_pass = false;}
17 
18  std::pair<bool, double> output(parameters_pass, weight);
19 
20  return output;
21 
22  }
output
BEGIN_PROLOG sequence::SlidingWindowTriggerPatternsOppositeWindows END_PROLOG simSlidingORM6O6 effSlidingORW output
Definition: triggersim_eastmodule_icarus.fcl:134
std::pair< bool, double > sbn::evwgh::FluxWeightCalc::PHSWCSVWeightCalc ( simb::MCFlux  flux,
std::vector< float >  rand 
)

Definition at line 10 of file PrimaryHadronSWCentralSplineVariationWeightCalc.cxx.

10  {
11 
12  //
13  // Largely built off the MiniBooNE code
14  // but this is intended to expand beyond it
15  //
16  // Edits from MiniBooNE code:
17  //
18  // JZ (6/2017) : Remove max weight, set to max_limit of a double
19  // JZ (6/2017) : Changed guards on c9 to be double point precision
20  //
21 
22  bool parameters_pass = true;
23 
24  double c1 = SWParam[0];
25  double c2 = SWParam[1];
26  double c3 = SWParam[2];
27  double c4 = SWParam[3];
28  double c5 = SWParam[4];
29  double c6 = SWParam[5];
30  double c7 = SWParam[6];
31  double c8 = SWParam[7];
32  double c9 = 1.0; // This isn't in the table but it is described in the text
33 
34  // Lay out the event kinimatics
35  double HadronMass = 0.13957010;
36 
37  // if(fabs(fprimaryHad) == 211) HadronMass = 0.13957010; //Charged Pion
38  // else{
39  // throw art::Exception(art::errors::StdException)
40  // << "sanford-wang is only configured for charged pions ";
41  // }
42 
43  TLorentzVector HadronVec;
44  double HadronPx = flux.ftpx;
45  double HadronPy = flux.ftpy;
46  double HadronPz = flux.ftpz;
47  double HadronE = sqrt(HadronPx*HadronPx +
48  HadronPy*HadronPy +
49  HadronPz*HadronPz +
50  HadronMass*HadronMass);
51  HadronVec.SetPxPyPzE(HadronPx,HadronPy,HadronPz,HadronE);
52 
53  ////////
54  //
55  // Based on MiniBooNE code to evaluate the theta value within below the maximum
56  // HARP theta coverage. This helps keep the splines well formed and constrains
57  // the uncertainties at very low neutrino energy
58  //
59  ////////
60 
61  double ThetaOfInterest;
62  if(HadronVec.Theta() > 0.195){
63  ThetaOfInterest = 0.195;
64  }
65  else{
66  ThetaOfInterest = HadronVec.Theta();
67  }
68 
69 
70  // Get Initial Proton Kinitmatics
71  // CURRENTLY GSimple flux files drop information about
72  // the initial state proton, but that this
73  TLorentzVector ProtonVec;
74  double ProtonMass = 0.9382720;
75  double ProtonPx = 0;
76  double ProtonPy = 0;
77  double ProtonPz = 8.89; //GeV
78  double ProtonE = sqrt(ProtonPx*ProtonPx +
79  ProtonPy*ProtonPy +
80  ProtonPz*ProtonPz +
81  ProtonMass*ProtonMass);
82  ProtonVec.SetPxPyPzE(ProtonPx,ProtonPy,ProtonPz,ProtonE);
83 
84  //Sanford-Wang Parameterization
85  // Eq 11 from PhysRevD.79.072002
86 
87  double CV = c1 * pow(HadronVec.P(), c2) *
88  (1. - HadronVec.P()/(ProtonVec.P() - c9)) *
89  exp(-1. * c3 * pow(HadronVec.P(), c4) / pow(ProtonVec.P(), c5)) *
90  exp(-1. * c6 * ThetaOfInterest *(HadronVec.P() - c7 * ProtonVec.P() * pow(cos(ThetaOfInterest), c8)));
91 
92  // Check taken from MiniBooNE code
93  if((HadronVec.P()) > ((ProtonVec.P()) - (c9))){
94  CV = 0;}
95 
96  //////
97  //
98  // Now that we have our central value based on the Sanford-Wang parameterization we
99  // need to create variations around this value. MiniBooNE did this by performing
100  // spline fits to the HARP data. This is a reproduction of that code:
101  //
102  ////
103 
104  double* HARPmomentumBins = HARPmomentumBounds.data(); // Convert std::vector to array
105  double* HARPthetaBins = HARPthetaBounds.data(); // Convert std::vector to array
106 
107  int Ntbins = int(HARPthetaBounds.size()) - 1;
108  int Npbins = int(HARPmomentumBounds.size()) - 1;
109 
110  //
111  // Using the HARP cross section and covariance matrices
112  // we will now create variations around the measured HARP
113  // meson production cross sections.
114  //
115 
116  // Important notes of the HARP data:
117  // The cross section matrix has 13 momentum bins and 6 theta bins
118  // The covariance matrix contains the correlated uncertainties across
119  // all 78 cross section measurements.
120  //
121  // The covariance matrix encodes the uncertainty on a given HARP
122  // ANALYSIS BIN instead of being a 3D matrix.
123  //
124  // A HARP analysis bin is defined as
125  // bin = momentum[theta[]] meaning that:
126  // analysis bin 0 is the zeroth theta and zeroth momentum bin
127  // but analysis bin 27 is the 3rd theta and 5th momentum bin
128  //
129  // The first thing to do is convert our cross section matrix into
130  // an std::vector for the analysis bins
131  //
132  std::vector< double > HARPCrossSectionAnalysisBins;
133  HARPCrossSectionAnalysisBins.resize(int(Ntbins*Npbins));
134 
135  int anaBin = 0;
136  for(int pbin = 0; pbin < Npbins; pbin++){
137  for(int tbin = 0; tbin < Ntbins; tbin++){
138  HARPCrossSectionAnalysisBins[anaBin] = HARPXSec[0][pbin][tbin];
139  anaBin++;
140  }
141  }
142 
143  //
144  // Now using this we can vary the cross section based on the HARP covariance matrix
145  // this will allow us to reweigh each cross section measurement based on
146  // the multigaussian smearing of this matrix.
147  //
148 
149  std::vector< double > smearedHARPCrossSectionAnalysisBins =
150  MultiGaussianSmearing(HARPCrossSectionAnalysisBins, FitCov, fIsDecomposed,rand);
151  HARPCrossSectionAnalysisBins.clear();
152 
153  //
154  // Check all the smeared cross sections, if any come out to be negative then
155  // we will not pass this given parameter set.
156  //
157 
158  for(int check = 0; check < int(smearedHARPCrossSectionAnalysisBins.size()); check++){
159  if(smearedHARPCrossSectionAnalysisBins[check] < 0){ parameters_pass = false;}
160  }
161 
162  //
163  // With a checked set of smeared cross sections we can convert our analysis bin
164  // std::vector back into a TMatrixD, though this is a bit annoying since TMatrices
165  // are an annoying format.
166  //
167  TMatrixD* smearedHARPXSec = new TMatrixD(Npbins,Ntbins);
168 
169  anaBin = 0;
170  for(int pbin = 0; pbin < Npbins; pbin++){
171  for(int tbin = 0; tbin < Ntbins; tbin++){
172  smearedHARPXSec[0][pbin][tbin] = smearedHARPCrossSectionAnalysisBins[anaBin];
173  anaBin++;
174  }
175  }
176 
177  //
178  // We want to now make a vector of histograms that are
179  // projections across the cross section matrix,
180  // this means we want to have vectors where one histogram
181  // is made for each bin on variable X where the histogram
182  // then contains bins for variable Y.
183  //
184  // This will allow us to do a 2D interpolation across the
185  // full parameter space using cubic spline fits
186  //
187 
188  std::vector< TH1F > MomentumBins;
189  MomentumBins.resize(Ntbins);
190  for(int bin = 0; bin < Ntbins; bin++){
191  MomentumBins[bin] = TH1F("HARPp",";;;", Npbins, HARPmomentumBins);
192  MomentumBins[bin].SetBit(kCanDelete);
193  }
194 
195  //Setup vectors of splines for fitting
196  std::vector< TSpline3 > SplinesVsMomentum;
197  SplinesVsMomentum.resize(Ntbins);
198 
199  for(int tbin = 0; tbin < Ntbins; tbin++){
200  for(int pbin = 0; pbin < Npbins; pbin++){
201 
202  //
203  // We want to create spline at fixed theta values, that way we can
204  // study the cross section at the given meson momentum
205  MomentumBins[tbin].SetBinContent(pbin+1, smearedHARPXSec[0][pbin][tbin]);
206  // This will give us the cross section in discreet slices of theta (tbin)
207  // but each histogram in the vector will be the cross section in bins of
208  // momentum that can then be cublic splined
209  //
210  }
211 
212  //
213  // in a given theta slice we can now spline over the cross section entries
214  // in bins of momentum
215  //
216  // Important note about Splines, MiniBooNE used constraints on the
217  // second derivative of the first and final knot point and required
218  // that they both equal to zero, this is not naturally the case in
219  // in TSpline3 but it is default in DCSPLC (the Fortran CERN library spline function)
220  // this helps to control the smoothness of the spline and minimizes the variation bin to bin
221  // For TSpine3 this is controlled by:
222  //
223  // b1 = constrain first knot 1st derivative
224  // b2 = constrain first knot 2nd derivative
225  // e1 = constrain final knot 1st derivative
226  // e2 = constrain final knot 2nd derivative
227  //
228  // the numbers that follow are the values you constrain those conditions to
229  //
230  SplinesVsMomentum[tbin] = TSpline3(&(MomentumBins[tbin]),"b2e2",0,0);
231 
232  }
233  MomentumBins.clear();
234  delete smearedHARPXSec;
235  //
236  // Now that we have the 1D Splines over momentum we want to know
237  // what the cross section is for the meson's momentum in bins of
238  // theta that way we can extract the cross section at the meson's
239  // theta.
240  //
241  TH1F* ThetaBins = new TH1F("HARPt",";;;", Ntbins, HARPthetaBins);
242  ThetaBins->SetBit(kCanDelete);
243  for(int tbin = 0; tbin < Ntbins; tbin++){
244  ThetaBins->SetBinContent(tbin+1, SplinesVsMomentum[tbin].Eval(HadronVec.P()));
245  }
246 
247  ////
248  // Now we can spline across the theta bins and we can extract the exact value of the
249  // cross section for the meson's theta value
250  //
251  // Options described above.
252  /////
253 
254  TSpline3* FinalSpline = new TSpline3(ThetaBins,"b2e2",0,0);
255 
256  double RW = FinalSpline->Eval(ThetaOfInterest);
257 
258  SplinesVsMomentum.clear();
259  delete ThetaBins;
260  delete FinalSpline;
261 
262  //
263  // These guards are inherited from MiniBooNE code
264  //
265  // These were defined here:
266  //
267  // cdcvs0.fnal.gov/cgi-bin/public-cvs/cvsweb-public.cgi/~checkout~/ ...
268  // miniboone/AnalysisFramework/MultisimMatrix/src/MultisimMatrix.inc
269  //
270  // This forces any negative spline fit to be 1
271 
272  ////////
273  // Possible Bug.
274  ////////
275  // This seems to be a feature in the MiniBooNE code
276  // It looks like the intension is to set this to zero
277  // but it is set to 1 before it is set to zero
278 
279  double weight = 1;
280 
281  if(RW < 0 || CV < 0){
282  weight = 1;
283  }
284  else if(fabs(CV) < 1.e-12){
285  weight = 1;
286  }
287  else{
288  weight *= RW/CV;
289  }
290 
291  if(weight < 0) weight = 0;
292  if(weight > 30) weight = 30;
293  if(!(std::isfinite(weight))){
294  std::cout << "SW+Splines : Failed to get a finite weight" << std::endl;
295  weight = 30;
296  }
297 
298  std::pair<bool, double> output(parameters_pass, weight);
299  return output;
300 
301  }// Done with the WeigthCalc function
sbn::evwgh::FluxWeightCalc::HARPmomentumBounds
std::vector< double > HARPmomentumBounds
Definition: FluxCalcPrep.h:76
icarus::ns::util::bin
constexpr details::BinObj< T > bin(T value)
Returns a wrapper to print the specified data in binary format.
geo::vect::isfinite
bool isfinite(Vector const &v)
Returns whether all components of the vector are finite.
Definition: geo_vectors_utils.h:2135
ThetaBins
BEGIN_PROLOG must be more than must be less than pi ThetaBins
Definition: GaisserParam.fcl:8
check
process_name tightIsolTest check
Definition: point_isolation_test.fcl:104
sbn::evwgh::FluxWeightCalc::SWParam
std::vector< double > SWParam
Definition: FluxCalcPrep.h:78
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
sbn::evwgh::FluxWeightCalc::HARPthetaBounds
std::vector< double > HARPthetaBounds
Definition: FluxCalcPrep.h:77
output
BEGIN_PROLOG sequence::SlidingWindowTriggerPatternsOppositeWindows END_PROLOG simSlidingORM6O6 effSlidingORW output
Definition: triggersim_eastmodule_icarus.fcl:134
sbn::evwgh::MultiGaussianSmearing
std::vector< std::vector< double > > MultiGaussianSmearing(std::vector< double > const &centralValue, std::vector< std::vector< double > > const &inputCovarianceMatrix, int n_multisims, CLHEP::RandGaussQ &GaussRandom)
Definition: SmearingUtils.cxx:10
e
do i e
Definition: updateVersionList.sh:68
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
std::pair< bool, double > sbn::evwgh::FluxWeightCalc::PHSWWeightCalc ( simb::MCFlux  flux,
std::vector< float >  rand 
)

Perform the MiniBooNE

Definition at line 9 of file PrimaryHadronSanfordWangWeightCalc.cxx.

9  {
10 
11  //
12  // Largely built off the MiniBooNE code
13  // but this is intended to expand beyond it
14  //
15  // Edits from MiniBooNE code:
16  //
17  // JZ (6/2017) : Remove max weight, set to max_limit of a double
18  // JZ (6/2017) : Changed guards on c9 to be double point precision
19  //
20 
21  // We need to guard against unphysical parameters
22  bool parameters_pass = true;
23 
24  // Get Neutrino Parent Kinimatics
25  double HadronMass = 0.4976;
26 
27 
28  // if(std::find(fprimaryHad.begin(), fprimaryHad.end(), 130) != fprimaryHad.end() ||
29  // std::find(fprimaryHad.begin(), fprimaryHad.end(), 310) != fprimaryHad.end() ||
30  // std::find(fprimaryHad.begin(), fprimaryHad.end(), 311) != fprimaryHad.end())
31  // {
32  // HadronMass = 0.4976; //Neutral Kaon
33  // }
34  // else{
35  // throw art::Exception(art::errors::StdException)
36  // << "Sanford-Wang is only configured for Netrual Kaons";
37  // }
38 
39  TLorentzVector HadronVec;
40  double HadronPx = flux.ftpx;
41  double HadronPy = flux.ftpy;
42  double HadronPz = flux.ftpz;
43  double HadronE = sqrt(HadronPx*HadronPx +
44  HadronPy*HadronPy +
45  HadronPz*HadronPz +
46  HadronMass*HadronMass);
47  HadronVec.SetPxPyPzE(HadronPx,HadronPy,HadronPz,HadronE);
48 
49  // Get Initial Proton Kinitmatics
50  // CURRENTLY GSimple flux files drop information about
51  // the initial state proton, but that this
52  TLorentzVector ProtonVec;
53  double ProtonMass = 0.9382720;
54  double ProtonPx = 0;
55  double ProtonPy = 0;
56  double ProtonPz = 8.89; //GeV
57  double ProtonE = sqrt(ProtonPx*ProtonPx +
58  ProtonPy*ProtonPy +
59  ProtonPz*ProtonPz +
60  ProtonMass*ProtonMass);
61  ProtonVec.SetPxPyPzE(ProtonPx,ProtonPy,ProtonPz,ProtonE);
62 
63  // Define the central value cross section
64  // Pull out the parameters for the Sanford-Wang Fit
65  double c1 = FitVal.at(0);
66  double c2 = FitVal.at(1);
67  double c3 = FitVal.at(2);
68  double c4 = FitVal.at(3);
69  double c5 = FitVal.at(4);
70  double c6 = FitVal.at(5);
71  double c7 = FitVal.at(6);
72  double c8 = FitVal.at(7);
73  double c9 = FitVal.at(8);
74 
75  //Sanford-Wang Parameterization
76  // Eq 11 from PhysRevD.79.072002
77 
78  double CV = c1 * pow(HadronVec.P(), c2) *
79  (1. - HadronVec.P()/(ProtonVec.P() - c9)) *
80  exp(-1. * c3 * pow(HadronVec.P(), c4) / pow(ProtonVec.P(), c5)) *
81  exp(-1. * c6 * HadronVec.Theta() *(HadronVec.P() - c7 * ProtonVec.P() * pow(cos(HadronVec.Theta()), c8)));
82 
83  // Check taken from MiniBooNE code
84  if((HadronVec.P()) > ((ProtonVec.P()) - (c9))){
85  CV = 0;
86  }
87 
88  // Define the variations around that cross section
89  // To do this we will call the a function from WeightCalc
90  // that will generate a set of smeared parameters based
91  // on the covariance matrix that we imported
92  std::vector< double > SWK0FitSmeared = MultiGaussianSmearing(FitVal, FitCov, rand);
93 
94  // Pull out the smeared parameters for the Sanford-Wang Fit
95  double smeared_c1 = SWK0FitSmeared.at(0);
96  double smeared_c2 = SWK0FitSmeared.at(1);
97  double smeared_c3 = SWK0FitSmeared.at(2);
98  double smeared_c4 = SWK0FitSmeared.at(3);
99  double smeared_c5 = SWK0FitSmeared.at(4);
100  double smeared_c6 = SWK0FitSmeared.at(5);
101  double smeared_c7 = SWK0FitSmeared.at(6);
102  double smeared_c8 = SWK0FitSmeared.at(7);
103  double smeared_c9 = SWK0FitSmeared.at(8);
104 
105  /// Perform the MiniBooNE
106  if(smeared_c1 < 0 ||
107  smeared_c3 < 0 ||
108  smeared_c6 < 0){
109  parameters_pass = false;
110  }
111 
112  double RW = smeared_c1 * pow(HadronVec.P(), smeared_c2) *
113  (1. - HadronVec.P()/(ProtonVec.P() - smeared_c9)) *
114  exp(-1. * smeared_c3 * pow(HadronVec.P(), smeared_c4) / pow(ProtonVec.P(), smeared_c5)) *
115  exp(-1. * smeared_c6 * HadronVec.Theta() *(HadronVec.P() - smeared_c7 * ProtonVec.P() * pow(cos(HadronVec.Theta()), smeared_c8)));
116 
117  // Check taken from MiniBooNE code
118  if((HadronVec.P()) > ((ProtonVec.P()) - (smeared_c9))){
119  RW = 0;
120  }
121 
122  double weight = 1;
123 
124  if(RW < 0 || CV < 0){//dont bother this; if this happens, the weight would be skipped...
125  weight = 1;
126  }
127  else if(fabs(CV) < 1.e-12){
128  weight = 1;
129  }
130  else{
131  weight *= RW/CV;
132  }
133 
134  if(weight < 0) weight = 0;
135  if(weight > 30) weight = 30;
136  if(!(std::isfinite(weight))){
137  std::cout << "SW : Failed to get a finite weight" << std::endl;
138  weight = 30;
139  }
140 
141  std::pair<bool, double> output(parameters_pass, weight);
142 
143  return output;
144 
145 
146  }
sbn::evwgh::FluxWeightCalc::FitVal
std::vector< double > FitVal
Definition: FluxCalcPrep.h:71
geo::vect::isfinite
bool isfinite(Vector const &v)
Returns whether all components of the vector are finite.
Definition: geo_vectors_utils.h:2135
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
output
BEGIN_PROLOG sequence::SlidingWindowTriggerPatternsOppositeWindows END_PROLOG simSlidingORM6O6 effSlidingORW output
Definition: triggersim_eastmodule_icarus.fcl:134
sbn::evwgh::MultiGaussianSmearing
std::vector< std::vector< double > > MultiGaussianSmearing(std::vector< double > const &centralValue, std::vector< std::vector< double > > const &inputCovarianceMatrix, int n_multisims, CLHEP::RandGaussQ &GaussRandom)
Definition: SmearingUtils.cxx:10
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
double sbn::evwgh::FluxWeightCalc::UnisimWeightCalc ( double  enu,
int  ptype,
int  ntype,
double  randomN,
bool  noNeg 
)

Definition at line 9 of file UnisimWeightCalc.cxx.

10  {//same copy from the FluxWeightCalc.cxx in ubsim/EventWeight/
11  //Keng Lin June 2021
12 
13  double weight = 1;
14 
15  int bin = int(enu/0.05); //convert energy (in GeV) into 50 MeV bin
16 
17 
18  // pseudocode:
19  // Scaled Reweighting = ScaleFactor * Reweighting + ( 1 - ScaleFactor) * Central Value
20  //
21  double scaled_pos = fScalePos*fRWpos[ptype][ntype][bin] +
22  (1-fScalePos)*fCV[ptype][ntype][bin];
23 
24  double scaled_neg = fScaleNeg*fRWneg[ptype][ntype][bin] +
25  (1-fScaleNeg)*fCV[ptype][ntype][bin];
26 
27  // pseudocode:
28  // Check value of Random Number for this universe such that:
29  // If random# > 0
30  // Weight = 1 + ( random# * [ { Scaled Reweighting[pos] / Central Value } - 1 ] )
31  // If random# < 0
32  // Weight = 1 - ( random# * [ { Scaled Reweighting[neg] / Central Value } - 1 ] )
33  //
34  // if there is only one systematic histogram offered sub this:
35  // If random# < 0
36  // Weight = 1 - ( random# * [ < 2 - { Scaled Reweighting[pos] / Central Value } > - 1 ] )
37  //
38 
39  if(randomN > 0){
40  double syst = randomN*((scaled_pos/fCV[ptype][ntype][bin])-1);
41  weight = 1 + (syst);
42 
43  if(scaled_pos == 0) weight = 1;
44  }
45 
46  else if(noNeg == true){
47  double syst = randomN*( (2 - (scaled_pos/fCV[ptype][ntype][bin])) - 1);
48  weight = 1 - (syst);
49 
50  if(scaled_pos == 0) weight = 1;
51 
52  }
53  else{
54  double syst = randomN*((scaled_neg/fCV[ptype][ntype][bin])-1);
55  weight = 1 - (syst);
56 
57  if(scaled_neg == 0) weight = 1;
58 
59  }
60 
61  if(fabs(fCV[ptype][ntype][bin]) < 1.e-12) weight = 1;
62 
63  if(weight < 0) weight = 1;
64  if(weight > 30) weight = 30;
65  if(!(std::isfinite(weight))){
66  std::cout << "UniSim " << fParameterSet.fName << " : Failed to get a finite weight" << std::endl;
67  weight = 30;}
68 
69  if( (ntype == 0 || ntype == 1) && ptype == 1) weight = 1;//nue/nuebar from pion
70  if( (ntype == 1 || ntype == 3) && ptype == 3) weight = 1;//nuebar/numubar from pi or charged kaon
71 
72  // std::cout<<weight<<" ("<<randomN<<") ";
73  return weight;
74 
75  }
icarus::ns::util::bin
constexpr details::BinObj< T > bin(T value)
Returns a wrapper to print the specified data in binary format.
sbn::evwgh::FluxWeightCalc::fRWpos
double fRWpos[4][4][200]
Definition: FluxCalcPrep.h:64
geo::vect::isfinite
bool isfinite(Vector const &v)
Returns whether all components of the vector are finite.
Definition: geo_vectors_utils.h:2135
ntype
then ntype
Definition: condor_hadd_sam.sh:202
sbn::evwgh::FluxWeightCalc::fRWneg
double fRWneg[4][4][200]
Definition: FluxCalcPrep.h:65
sbn::evwgh::EventWeightParameterSet::fName
std::string fName
Name of the parameter set.
Definition: EventWeightParameterSet.h:145
sbn::evwgh::FluxWeightCalc::fCV
double fCV[4][4][200]
Definition: FluxCalcPrep.h:63
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10

Member Data Documentation

std::string sbn::evwgh::FluxWeightCalc::CalcType
private

Definition at line 56 of file FluxCalcPrep.h.

double sbn::evwgh::FluxWeightCalc::fCV[4][4][200]
private

Definition at line 63 of file FluxCalcPrep.h.

std::string sbn::evwgh::FluxWeightCalc::fGeneratorModuleLabel
private

Definition at line 55 of file FluxCalcPrep.h.

bool sbn::evwgh::FluxWeightCalc::fIsDecomposed {false}
private

Definition at line 79 of file FluxCalcPrep.h.

TMatrixD* sbn::evwgh::FluxWeightCalc::FitCov {nullptr}
private

Definition at line 72 of file FluxCalcPrep.h.

std::vector<double> sbn::evwgh::FluxWeightCalc::FitVal {}
private

Definition at line 71 of file FluxCalcPrep.h.

std::vector<int> sbn::evwgh::FluxWeightCalc::fprimaryHad = {0}
private

Definition at line 69 of file FluxCalcPrep.h.

double sbn::evwgh::FluxWeightCalc::fRWneg[4][4][200]
private

Definition at line 65 of file FluxCalcPrep.h.

double sbn::evwgh::FluxWeightCalc::fRWpos[4][4][200]
private

Definition at line 64 of file FluxCalcPrep.h.

double sbn::evwgh::FluxWeightCalc::fScaleNeg = 1
private

Definition at line 59 of file FluxCalcPrep.h.

double sbn::evwgh::FluxWeightCalc::fScalePos {}
private

Definition at line 58 of file FluxCalcPrep.h.

std::vector<double> sbn::evwgh::FluxWeightCalc::HARPmomentumBounds {}
private

Definition at line 76 of file FluxCalcPrep.h.

std::vector<double> sbn::evwgh::FluxWeightCalc::HARPthetaBounds {}
private

Definition at line 77 of file FluxCalcPrep.h.

TMatrixD* sbn::evwgh::FluxWeightCalc::HARPXSec {nullptr}
private

Definition at line 75 of file FluxCalcPrep.h.

bool sbn::evwgh::FluxWeightCalc::PosOnly {false}
private

Definition at line 66 of file FluxCalcPrep.h.

std::vector<double> sbn::evwgh::FluxWeightCalc::SWParam {}
private

Definition at line 78 of file FluxCalcPrep.h.

int sbn::evwgh::FluxWeightCalc::wc = 0
private

Definition at line 82 of file FluxCalcPrep.h.

int sbn::evwgh::FluxWeightCalc::wc0 = 0
private

Definition at line 84 of file FluxCalcPrep.h.

int sbn::evwgh::FluxWeightCalc::wc1 = 0
private

Definition at line 85 of file FluxCalcPrep.h.

int sbn::evwgh::FluxWeightCalc::wc30 = 0
private

Definition at line 86 of file FluxCalcPrep.h.

int sbn::evwgh::FluxWeightCalc::wcn = 0
private

Definition at line 83 of file FluxCalcPrep.h.

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