00001 #ifndef neugenweightcalculator_cxx
00002 #define neugenweightcalculator_cxx
00003 #include "MessageService/MsgService.h"
00004 #include "MCReweight/NeugenWeightCalculator.h"
00005 #include "NeugenInterface/neugen_wrapper.h"
00006 #include "NeugenInterface/neugen_config.h"
00007 #include "NeugenInterface/interaction.h"
00008 #include "NeugenInterface/process.h"
00009 #include "NeugenInterface/kinematic_variable.h"
00010 #include "NeugenInterface/init_state.h"
00011 #include "NeugenInterface/final_state.h"
00012 #include "NeugenInterface/flavor.h"
00013 #include "TLorentzVector.h"
00014 #include <cassert>
00015
00016 CVSID("$Id: NeugenWeightCalculator.cxx,v 1.22 2008/02/22 21:48:43 boehm Exp $");
00017
00018
00019 NeugenWeightCalculator::NeugenWeightCalculator(Registry *config)
00020 : WeightCalculator(config)
00021 {
00022
00023 fWCname = "NeugenWeightCalculator";
00024
00025 fStandardConfigChanged = false;
00026 fReweightConfigChanged = false;
00027
00028 fStdNeuConfig = new neugen_config("stdNeuConfig");
00029 fStdNeuConfig->set_best_parameters();
00030 if(config) SetStandardConfig(config);
00031
00032 fRwtNeuConfig = new neugen_config("rwtNeuConfig");
00033 fRwtNeuConfig->set_best_parameters();
00034 if(config) SetReweightConfig(config);
00035
00036 fNuWrap = new neugen_wrapper(fStdNeuConfig);
00037
00038 fIsOptimized = false;
00039 fUpdateConfig = false;
00040 }
00041
00042
00043 NeugenWeightCalculator::~NeugenWeightCalculator()
00044 {
00045 delete fStdNeuConfig;
00046 delete fRwtNeuConfig;
00047 delete fNuWrap;
00048 }
00049
00050
00051 void NeugenWeightCalculator::ReweightConfigReset()
00052 {
00053 fRwtNeuConfig->set_best_parameters();
00054 if(fStandardConfig) SetReweightConfig(fStandardConfig);
00055 }
00056
00057
00058
00059 void NeugenWeightCalculator::SetOptimization(bool in)
00060 {
00061 fIsOptimized = in;
00062 }
00063
00064
00065 void NeugenWeightCalculator::Config()
00066 {
00067
00068 double tempd = 0;
00069 int tempi = 0;
00070 const char *temps =0;
00071
00072 if(fStandardConfigChanged) {
00073 bool setname = false;
00074 bool setno = false;
00075 if(fStandardConfig->Get("neugen:config_no",tempi)){
00076 fStdNeuConfig->set_config_no(tempi);
00077 setno = true;
00078 }
00079 if(fStandardConfig->Get("neugen:config_name",temps)){
00080 fStdNeuConfig->set_config_name(temps);
00081 setname=true;
00082 }
00083 if(setname&&!setno){
00084 std::cout<<"You tried to set a new config name in NeugenWeightCalculator::Config()"
00085 <<" but you did not set a config no. Dire consequences will most likely result"
00086 <<std::endl;
00087 assert(0);
00088 }
00089
00090
00091
00092
00093
00094
00095 if(setname&&setno){
00096 fStdNeuConfig->fill_named_configuration();
00097 }
00098
00099 if(fStandardConfig->Get("neugen:pdf_group",tempi))
00100 fStdNeuConfig->set_pdf_group(tempi);
00101 if(fStandardConfig->Get("neugen:pdf_set",tempi))
00102 fStdNeuConfig->set_pdf_set(tempi);
00103 if(fStandardConfig->Get("neugen:ma_qe",tempd))
00104 fStdNeuConfig->set_ma_qe(tempd);
00105 if(fStandardConfig->Get("neugen:ma_res",tempd))
00106 fStdNeuConfig->set_ma_res(tempd);
00107 if(fStandardConfig->Get("neugen:ma_coh",tempd))
00108 fStdNeuConfig->set_ma_coh(tempd);
00109 if(fStandardConfig->Get("neugen:qel_fa0",tempd))
00110 fStdNeuConfig->set_qel_fa0(tempd);
00111 if(fStandardConfig->Get("neugen:qel_eta",tempd))
00112 fStdNeuConfig->set_qel_eta(tempd);
00113 if(fStandardConfig->Get("neugen:res_omega",tempd))
00114 fStdNeuConfig->set_res_omega(tempd);
00115 if(fStandardConfig->Get("neugen:res_z",tempd))
00116 fStdNeuConfig->set_res_z(tempd);
00117 if(fStandardConfig->Get("neugen:coh_r0",tempd))
00118 fStdNeuConfig->set_coh_r0(tempd);
00119 if(fStandardConfig->Get("neugen:coh_rei",tempd))
00120 fStdNeuConfig->set_coh_rei(tempd);
00121 if(fStandardConfig->Get("neugen:kno_a1",tempd))
00122 fStdNeuConfig->set_kno_a1(tempd);
00123 if(fStandardConfig->Get("neugen:kno_a2",tempd))
00124 fStdNeuConfig->set_kno_a2(tempd);
00125 if(fStandardConfig->Get("neugen:kno_a3",tempd))
00126 fStdNeuConfig->set_kno_a3(tempd);
00127 if(fStandardConfig->Get("neugen:kno_a4",tempd))
00128 fStdNeuConfig->set_kno_a4(tempd);
00129 if(fStandardConfig->Get("neugen:kno_b1",tempd))
00130 fStdNeuConfig->set_kno_b1(tempd);
00131 if(fStandardConfig->Get("neugen:kno_b2",tempd))
00132 fStdNeuConfig->set_kno_b2(tempd);
00133 if(fStandardConfig->Get("neugen:kno_b3",tempd))
00134 fStdNeuConfig->set_kno_b3(tempd);
00135 if(fStandardConfig->Get("neugen:kno_b4",tempd))
00136 fStdNeuConfig->set_kno_b4(tempd);
00137 if(fStandardConfig->Get("neugen:kno_c1",tempd))
00138 fStdNeuConfig->set_kno_c1(tempd);
00139 if(fStandardConfig->Get("neugen:kno_c2",tempd))
00140 fStdNeuConfig->set_kno_c2(tempd);
00141 if(fStandardConfig->Get("neugen:kno_c3",tempd))
00142 fStdNeuConfig->set_kno_c3(tempd);
00143 if(fStandardConfig->Get("neugen:kno_c4",tempd))
00144 fStdNeuConfig->set_kno_c4(tempd);
00145
00146
00147 if(fStandardConfig->Get("neugen:kno_r112",tempd))
00148 fStdNeuConfig->set_dis_res(1,2,e_vp,tempd);
00149 if(fStandardConfig->Get("neugen:kno_r122",tempd))
00150 fStdNeuConfig->set_dis_res(1,2,e_vn,tempd);
00151 if(fStandardConfig->Get("neugen:kno_r132",tempd))
00152 fStdNeuConfig->set_dis_res(1,2,e_vbp,tempd);
00153 if(fStandardConfig->Get("neugen:kno_r142",tempd))
00154 fStdNeuConfig->set_dis_res(1,2,e_vbn,tempd);
00155 if(fStandardConfig->Get("neugen:kno_r113",tempd))
00156 fStdNeuConfig->set_dis_res(1,3,e_vp,tempd);
00157 if(fStandardConfig->Get("neugen:kno_r123",tempd))
00158 fStdNeuConfig->set_dis_res(1,3,e_vn,tempd);
00159 if(fStandardConfig->Get("neugen:kno_r133",tempd))
00160 fStdNeuConfig->set_dis_res(1,3,e_vbp,tempd);
00161 if(fStandardConfig->Get("neugen:kno_r143",tempd))
00162 fStdNeuConfig->set_dis_res(1,3,e_vbn,tempd);
00163
00164 if(fStandardConfig->Get("neugen:kno_r212",tempd))
00165 fStdNeuConfig->set_dis_res(2,2,e_vp,tempd);
00166 if(fStandardConfig->Get("neugen:kno_r222",tempd))
00167 fStdNeuConfig->set_dis_res(2,2,e_vn,tempd);
00168 if(fStandardConfig->Get("neugen:kno_r232",tempd))
00169 fStdNeuConfig->set_dis_res(2,2,e_vbp,tempd);
00170 if(fStandardConfig->Get("neugen:kno_r242",tempd))
00171 fStdNeuConfig->set_dis_res(2,2,e_vbn,tempd);
00172 if(fStandardConfig->Get("neugen:kno_r213",tempd))
00173 fStdNeuConfig->set_dis_res(2,3,e_vp,tempd);
00174 if(fStandardConfig->Get("neugen:kno_r223",tempd))
00175 fStdNeuConfig->set_dis_res(2,3,e_vn,tempd);
00176 if(fStandardConfig->Get("neugen:kno_r233",tempd))
00177 fStdNeuConfig->set_dis_res(2,3,e_vbp,tempd);
00178 if(fStandardConfig->Get("neugen:kno_r243",tempd))
00179 fStdNeuConfig->set_dis_res(2,3,e_vbn,tempd);
00180
00181 if(fStandardConfig->Get("neugen:wcutd",tempd))
00182 fStdNeuConfig->set_wcutd(tempd);
00183 if(fStandardConfig->Get("neugen:wcutr",tempd))
00184 fStdNeuConfig->set_wcutr(tempd);
00185 if(fStandardConfig->Get("neugen:nres",tempi))
00186 fStdNeuConfig->set_nres(tempi);
00187 if(fStandardConfig->Get("neugen:dcf",tempd))
00188 fStdNeuConfig->set_dcf(tempd);
00189
00190 fStandardConfigChanged = false;
00191 }
00192
00193
00194 if(fReweightConfigChanged) {
00195 fUpdateConfig = true;
00196
00197
00198 Float_t sc = 0; tempi = 0;
00199 if(fReweightConfig->Get("neugen:use_scale_factors",tempi)){
00200 if(tempi==1) sc = 1.0;
00201 }
00202
00203
00204
00205
00206 bool gotname=false;
00207 bool gotno=false;
00208 if(fReweightConfig->Get("neugen:config_name",temps)){
00209 fRwtNeuConfig->set_config_name(temps);
00210 gotname=true;
00211 }
00212
00213 if(fReweightConfig->Get("neugen:config_no",tempi)){
00214 fRwtNeuConfig->set_config_no(tempi);
00215 gotno=true;
00216 }
00217
00218 if(!gotname||!gotno){
00219 fRwtNeuConfig->set_config_name(fStdNeuConfig->config_name());
00220 fRwtNeuConfig->set_config_no(fStdNeuConfig->config_no());
00221 }
00222
00223
00224
00225 fRwtNeuConfig->fill_named_configuration();
00226
00227
00228
00229
00230 if(fReweightConfig->Get("neugen:pdf_group",tempi))
00231 fRwtNeuConfig->set_pdf_group(tempi);
00232 if(fReweightConfig->Get("neugen:pdf_set",tempi))
00233 fRwtNeuConfig->set_pdf_set(tempi);
00234
00235 if(fReweightConfig->Get("neugen:ma_qe",tempd))
00236 fRwtNeuConfig->set_ma_qe(tempd*(1+sc*(fRwtNeuConfig->ma_qe()-1)));
00237 if(fReweightConfig->Get("neugen:ma_res",tempd))
00238 fRwtNeuConfig->set_ma_res(tempd*(1+sc*(fRwtNeuConfig->ma_res()-1)));
00239 if(fReweightConfig->Get("neugen:ma_coh",tempd))
00240 fRwtNeuConfig->set_ma_coh(tempd*(1+sc*(fRwtNeuConfig->ma_coh()-1)));
00241 if(fReweightConfig->Get("neugen:qel_fa0",tempd))
00242 fRwtNeuConfig->set_qel_fa0(tempd*(1+sc*(fRwtNeuConfig->qel_fa0()-1)));
00243 if(fReweightConfig->Get("neugen:qel_eta",tempd))
00244 fRwtNeuConfig->set_qel_eta(tempd*(1+sc*(fRwtNeuConfig->qel_eta()-1)));
00245 if(fReweightConfig->Get("neugen:res_omega",tempd))
00246 fRwtNeuConfig->set_res_omega(tempd*(1+sc*(fRwtNeuConfig->res_omega()-1)));
00247 if(fReweightConfig->Get("neugen:res_z",tempd))
00248 fRwtNeuConfig->set_res_z(tempd*(1+sc*(fRwtNeuConfig->res_z()-1)));
00249 if(fReweightConfig->Get("neugen:coh_r0",tempd))
00250 fRwtNeuConfig->set_coh_r0(tempd*(1+sc*(fRwtNeuConfig->coh_r0()-1)));
00251 if(fReweightConfig->Get("neugen:coh_rei",tempd))
00252 fRwtNeuConfig->set_coh_rei(tempd*(1+sc*(fRwtNeuConfig->coh_rei()-1)));
00253 if(fReweightConfig->Get("neugen:kno_a1",tempd))
00254 fRwtNeuConfig->set_kno_a1(tempd*(1+sc*(fRwtNeuConfig->kno_a1()-1)));
00255 if(fReweightConfig->Get("neugen:kno_a2",tempd))
00256 fRwtNeuConfig->set_kno_a2(tempd*(1+sc*(fRwtNeuConfig->kno_a2()-1)));
00257 if(fReweightConfig->Get("neugen:kno_a3",tempd))
00258 fRwtNeuConfig->set_kno_a3(tempd*(1+sc*(fRwtNeuConfig->kno_a3()-1)));
00259 if(fReweightConfig->Get("neugen:kno_a4",tempd))
00260 fRwtNeuConfig->set_kno_a4(tempd*(1+sc*(fRwtNeuConfig->kno_a4()-1)));
00261 if(fReweightConfig->Get("neugen:kno_b1",tempd))
00262 fRwtNeuConfig->set_kno_b1(tempd*(1+sc*(fRwtNeuConfig->kno_b1()-1)));
00263 if(fReweightConfig->Get("neugen:kno_b2",tempd))
00264 fRwtNeuConfig->set_kno_b2(tempd*(1+sc*(fRwtNeuConfig->kno_b2()-1)));
00265 if(fReweightConfig->Get("neugen:kno_b3",tempd))
00266 fRwtNeuConfig->set_kno_b3(tempd*(1+sc*(fRwtNeuConfig->kno_b3()-1)));
00267 if(fReweightConfig->Get("neugen:kno_b4",tempd))
00268 fRwtNeuConfig->set_kno_b4(tempd*(1+sc*(fRwtNeuConfig->kno_b4()-1)));
00269 if(fReweightConfig->Get("neugen:kno_c1",tempd))
00270 fRwtNeuConfig->set_kno_c1(tempd*(1+sc*(fRwtNeuConfig->kno_c1()-1)));
00271 if(fReweightConfig->Get("neugen:kno_c2",tempd))
00272 fRwtNeuConfig->set_kno_c2(tempd*(1+sc*(fRwtNeuConfig->kno_c2()-1)));
00273 if(fReweightConfig->Get("neugen:kno_c3",tempd))
00274 fRwtNeuConfig->set_kno_c3(tempd*(1+sc*(fRwtNeuConfig->kno_c3()-1)));
00275 if(fReweightConfig->Get("neugen:kno_c4",tempd))
00276 fRwtNeuConfig->set_kno_c4(tempd*(1+sc*(fRwtNeuConfig->kno_c4()-1)));
00277
00278
00279
00280 if(fReweightConfig->Get("neugen:kno_r112",tempd))
00281 fRwtNeuConfig->set_dis_res(1,2,e_vp,tempd*(1+sc*(fRwtNeuConfig->dis_res(1,2,e_vp)-1)));
00282 if(fReweightConfig->Get("neugen:kno_r122",tempd))
00283 fRwtNeuConfig->set_dis_res(1,2,e_vn,tempd*(1+sc*(fRwtNeuConfig->dis_res(1,2,e_vn)-1)));
00284 if(fReweightConfig->Get("neugen:kno_r132",tempd))
00285 fRwtNeuConfig->set_dis_res(1,2,e_vbp,tempd*(1+sc*(fRwtNeuConfig->dis_res(1,2,e_vbp)-1)));
00286 if(fReweightConfig->Get("neugen:kno_r142",tempd))
00287 fRwtNeuConfig->set_dis_res(1,2,e_vbn,tempd*(1+sc*(fRwtNeuConfig->dis_res(1,2,e_vbn)-1)));
00288 if(fReweightConfig->Get("neugen:kno_r113",tempd))
00289 fRwtNeuConfig->set_dis_res(1,3,e_vp,tempd*(1+sc*(fRwtNeuConfig->dis_res(1,3,e_vp)-1)));
00290 if(fReweightConfig->Get("neugen:kno_r123",tempd))
00291 fRwtNeuConfig->set_dis_res(1,3,e_vn,tempd*(1+sc*(fRwtNeuConfig->dis_res(1,3,e_vn)-1)));
00292 if(fReweightConfig->Get("neugen:kno_r133",tempd))
00293 fRwtNeuConfig->set_dis_res(1,3,e_vbp,tempd*(1+sc*(fRwtNeuConfig->dis_res(1,3,e_vbp)-1)));
00294 if(fReweightConfig->Get("neugen:kno_r143",tempd))
00295 fRwtNeuConfig->set_dis_res(1,3,e_vbn,tempd*(1+sc*(fRwtNeuConfig->dis_res(1,3,e_vbn)-1)));
00296
00297 if(fReweightConfig->Get("neugen:kno_r212",tempd))
00298 fRwtNeuConfig->set_dis_res(2,2,e_vp,tempd*(1+sc*(fRwtNeuConfig->dis_res(2,2,e_vp)-1)));
00299 if(fReweightConfig->Get("neugen:kno_r222",tempd))
00300 fRwtNeuConfig->set_dis_res(2,2,e_vn,tempd*(1+sc*(fRwtNeuConfig->dis_res(2,2,e_vn)-1)));
00301 if(fReweightConfig->Get("neugen:kno_r232",tempd))
00302 fRwtNeuConfig->set_dis_res(2,2,e_vbp,tempd*(1+sc*(fRwtNeuConfig->dis_res(2,2,e_vbp)-1)));
00303 if(fReweightConfig->Get("neugen:kno_r242",tempd))
00304 fRwtNeuConfig->set_dis_res(2,2,e_vbn,tempd*(1+sc*(fRwtNeuConfig->dis_res(2,2,e_vbn)-1)));
00305 if(fReweightConfig->Get("neugen:kno_r213",tempd))
00306 fRwtNeuConfig->set_dis_res(2,3,e_vp,tempd*(1+sc*(fRwtNeuConfig->dis_res(2,3,e_vp)-1)));
00307 if(fReweightConfig->Get("neugen:kno_r223",tempd))
00308 fRwtNeuConfig->set_dis_res(2,3,e_vn,tempd*(1+sc*(fRwtNeuConfig->dis_res(2,3,e_vn)-1)));
00309 if(fReweightConfig->Get("neugen:kno_r233",tempd))
00310 fRwtNeuConfig->set_dis_res(2,3,e_vbp,tempd*(1+sc*(fRwtNeuConfig->dis_res(2,3,e_vbp)-1)));
00311 if(fReweightConfig->Get("neugen:kno_r243",tempd))
00312 fRwtNeuConfig->set_dis_res(2,3,e_vbn,tempd*(1+sc*(fRwtNeuConfig->dis_res(2,3,e_vbn)-1)));
00313
00314
00315 if(fReweightConfig->Get("neugen:scale_kno_all",tempd)){
00316 for(int c=1;c<=2;c++){
00317 for(int m=2;m<=3;m++){
00318 fRwtNeuConfig->set_dis_res(c,m,e_vp,
00319 fRwtNeuConfig->dis_res(c,m,e_vp)*tempd);
00320 fRwtNeuConfig->set_dis_res(c,m,e_vn,
00321 fRwtNeuConfig->dis_res(c,m,e_vn)*tempd);
00322 fRwtNeuConfig->set_dis_res(c,m,e_vbp,
00323 fRwtNeuConfig->dis_res(c,m,e_vbp)*tempd);
00324 fRwtNeuConfig->set_dis_res(c,m,e_vbn,
00325 fRwtNeuConfig->dis_res(c,m,e_vbn)*tempd);
00326 }
00327 }
00328 }
00329 if(fReweightConfig->Get("neugen:scale_kno_cc",tempd)){
00330 for(int m=2;m<=3;m++){
00331 fRwtNeuConfig->set_dis_res(1,m,e_vp,
00332 fRwtNeuConfig->dis_res(1,m,e_vp)*tempd);
00333 fRwtNeuConfig->set_dis_res(1,m,e_vn,
00334 fRwtNeuConfig->dis_res(1,m,e_vn)*tempd);
00335 fRwtNeuConfig->set_dis_res(1,m,e_vbp,
00336 fRwtNeuConfig->dis_res(1,m,e_vbp)*tempd);
00337 fRwtNeuConfig->set_dis_res(1,m,e_vbn,
00338 fRwtNeuConfig->dis_res(1,m,e_vbn)*tempd);
00339 }
00340 }
00341 if(fReweightConfig->Get("neugen:scale_kno_nc",tempd)){
00342 for(int m=2;m<=3;m++){
00343 fRwtNeuConfig->set_dis_res(2,m,e_vp,
00344 fRwtNeuConfig->dis_res(2,m,e_vp)*tempd);
00345 fRwtNeuConfig->set_dis_res(2,m,e_vn,
00346 fRwtNeuConfig->dis_res(2,m,e_vn)*tempd);
00347 fRwtNeuConfig->set_dis_res(2,m,e_vbp,
00348 fRwtNeuConfig->dis_res(2,m,e_vbp)*tempd);
00349 fRwtNeuConfig->set_dis_res(2,m,e_vbn,
00350 fRwtNeuConfig->dis_res(2,m,e_vbn)*tempd);
00351 }
00352 }
00353 if(fReweightConfig->Get("neugen:wcutd",tempd))
00354 fRwtNeuConfig->set_wcutd(tempd*(1+sc*(fRwtNeuConfig->wcutd()-1)));
00355 if(fReweightConfig->Get("neugen:wcutr",tempd))
00356 fRwtNeuConfig->set_wcutr(tempd*(1+sc*(fRwtNeuConfig->wcutr()-1)));
00357 if(fReweightConfig->Get("neugen:nres",tempi))
00358 fRwtNeuConfig->set_nres(tempi);
00359 if(fReweightConfig->Get("neugen:dcf",tempd))
00360 fRwtNeuConfig->set_dcf(tempd*(1+sc*(fRwtNeuConfig->dcf()-1)));
00361
00362
00363
00364 fReweightConfigChanged = false;
00365 }
00366 }
00367
00368
00369 double NeugenWeightCalculator::GetWeight(Registry *event)
00370 {
00371 if(!event) return 1;
00372
00373 double nu_en = 0;
00374 double nu_px = 0;
00375 double nu_py = 0;
00376 double nu_pz = 0;
00377 if(!event->Get("event:nu_en",nu_en)) return 1;
00378 if(!event->Get("event:nu_px",nu_px)) return 1;
00379 if(!event->Get("event:nu_py",nu_py)) return 1;
00380 if(!event->Get("event:nu_pz",nu_pz)) return 1;
00381
00382
00383 double tar_en = 0;
00384 double tar_px = 0;
00385 double tar_py = 0;
00386 double tar_pz = 0;
00387 if(!event->Get("event:tar_en",tar_en)) return 1;
00388 if(!event->Get("event:tar_px",tar_px)) return 1;
00389 if(!event->Get("event:tar_py",tar_py)) return 1;
00390 if(!event->Get("event:tar_pz",tar_pz)) return 1;
00391
00392
00393 double kval[5] = {0,0,0,0,0};
00394 int kid1 = -1;
00395 int kid2 = -1;
00396 if(!event->Get("event:q2",kval[1])) kval[1] = 999999;
00397 if(!event->Get("event:w2",kval[2])) kval[2] = 999999;
00398 if(!event->Get("event:x",kval[3])) kval[3] = 999999;
00399 if(!event->Get("event:y",kval[4])) kval[4] = 999999;
00400
00401
00402 int iaction = 0;
00403 int inu = 0;
00404 int ires = 0;
00405 int initial_state = 0;
00406 int nucleus = 0;
00407 int had_fs = 0;
00408 if(!event->Get("event:iaction",iaction)) return 1;
00409 if(!event->Get("event:inu",inu)) return 1;
00410 if(!event->Get("event:iresonance",ires)) {
00411 if(!event->Get("event:process",ires)) return 1;
00412 }
00413 if(!event->Get("event:initial_state",initial_state)) return 1;
00414 if(!event->Get("event:nucleus",nucleus)) return 1;
00415 if(!event->Get("event:had_fs",had_fs)) return 1;
00416
00417
00418 int flavor = 4;
00419 if(abs(inu)==12) flavor = 1;
00420 else if(abs(inu)==14) flavor = 2;
00421 else if(abs(inu)==16) flavor = 3;
00422 else if(abs(inu)==1 || abs(inu)==2 || abs(inu)==3) flavor = abs(inu);
00423 else {
00424 MSG("NeugenWC",Msg::kError) << "Invalid event:inu value " << inu << endl;
00425 return 1;
00426 }
00427
00428 int ccnc = iaction;
00429 if(ccnc==0) ccnc = 2;
00430 if(!(ccnc==1||ccnc==2)) {
00431 MSG("NeugenWC",Msg::kError) << "Invalid event:iaction value "
00432 << iaction << endl;
00433 return 1;
00434 }
00435
00436 int process = 0;
00437 if(ires>=1001&&ires<=1004) process = ires - 1000;
00438 else if(ires>=1&&ires<=4) process = ires;
00439 else {
00440 if(ires==1005) MSG("NeugenWC",Msg::kInfo)
00441 << "Event:iresonance value is 1005 "
00442 << "- no reweighting performed" << endl;
00443 else MSG("NeugenWC",Msg::kError)
00444 << "Invalid event:iresonance (or event:process) value "
00445 << ires << endl;
00446 return 1;
00447 }
00448
00449
00450 if(process==1) {
00451 kid1=1;
00452 kid2=0;
00453 if(kval[1]==999999) {
00454 MSG("NeugenWC",Msg::kWarning)
00455 << "Current event is QEL: q^2 kinematic variable not set!"
00456 << " Returning weight of 1." << endl;
00457 return 1;
00458 }
00459 }
00460 else if(process==2) {
00461 kid1=1;
00462 kid2=2;
00463 if(kval[1]==999999) {
00464 MSG("NeugenWC",Msg::kWarning)
00465 << "Current event is RES: q^2 kinematic variable not set!"
00466 << " Returning weight of 1." << endl;
00467 return 1;
00468 }
00469 if(kval[2]==999999) {
00470 MSG("NeugenWC",Msg::kWarning)
00471 << "Current event is RES: w^2 kinematic variable not set!"
00472 << " Returning weight of 1." << endl;
00473 return 1;
00474 }
00475
00476 if(kval[2]>0) kval[2]=sqrt(kval[2]);
00477 else {
00478 MSG("NeugenWC",Msg::kWarning)
00479 << "w^2 kinematic variable is negative!"
00480 << " Returning weight of 1." << endl;
00481 return 1;
00482 }
00483 }
00484 else if(process==3) {
00485 kid1=3;
00486 kid2=4;
00487 if(kval[3]==999999) {
00488 MSG("NeugenWC",Msg::kWarning)
00489 << "Current event is DIS: x kinematic variable not set!"
00490 << " Returning weight of 1." << endl;
00491 return 1;
00492 }
00493 if(kval[4]==999999) {
00494 MSG("NeugenWC",Msg::kWarning)
00495 << "Current event is DIS: y kinematic variable not set!"
00496 << " Returning weight of 1." << endl;
00497 return 1;
00498 }
00499 }
00500 else if(process==4) {
00501
00502 return 1;
00503 }
00504
00505
00506 TLorentzVector *nu_lv = new TLorentzVector(nu_px,nu_py,nu_pz,nu_en);
00507 TLorentzVector *tar_lv = new TLorentzVector(tar_px,tar_py,tar_pz,tar_en);
00508 interaction *iact = new interaction(flavor_enum(flavor),
00509 nucleus_enum(nucleus),
00510 ccnc_enum(ccnc),
00511 init_state_enum(initial_state));
00512 final_state *fs = new final_state(had_fs);
00513
00514
00515 double weight = 0;
00516 weight = fNuWrap->reweight(nu_lv,tar_lv,
00517 kinematic_variable_enum(kid1),float(kval[kid1]),
00518 kinematic_variable_enum(kid2),float(kval[kid2]),
00519 iact,process_enum(process),fs,
00520 fRwtNeuConfig);
00521
00522 delete nu_lv;
00523 delete tar_lv;
00524 delete iact;
00525 delete fs;
00526 return weight;
00527
00528 }
00529
00530
00531 double NeugenWeightCalculator::GetWeight(MCEventInfo *event,NuParent *parent)
00532 {
00533
00534 if(!event) return 1;
00535 if(!parent) {
00536 MSG("NeugenWC",Msg::kDebug) << "No NuParent object" << endl;
00537 }
00538
00539
00540 double nu_en = event->nuE;
00541 double nu_px = event->nuPx;
00542 double nu_py = event->nuPy;
00543 double nu_pz = event->nuPz;
00544 if(nu_en==0) return 1;
00545
00546
00547 double tar_en = event->tarE;
00548 double tar_px = event->tarPx;
00549 double tar_py = event->tarPy;
00550 double tar_pz = event->tarPz;
00551 if(tar_en==0) return 1;
00552
00553
00554 double kval[5] = {0,0,0,0,0};
00555 int kid1 = -1;
00556 int kid2 = -1;
00557 kval[1] = event->q2;
00558 kval[2] = event->w2;
00559 kval[3] = event->x;
00560 kval[4] = event->y;
00561
00562
00563 int iaction = event->iaction;
00564 int inu = event->inu;
00565 int ires = event->iresonance;
00566 int initial_state = event->initial_state;
00567 int nucleus = event->nucleus;
00568 int had_fs = event->had_fs;
00569
00570
00571 int flavor = 4;
00572 if(abs(inu)==12) flavor = 1;
00573 else if(abs(inu)==14) flavor = 2;
00574 else if(abs(inu)==16) flavor = 3;
00575 else if(abs(inu)==1 || abs(inu)==2 || abs(inu)==3) flavor = abs(inu);
00576 else {
00577 MSG("NeugenWC",Msg::kError) << "Invalid event:inu value " << inu << endl;
00578 return 1;
00579 }
00580
00581 int ccnc = iaction;
00582 if(ccnc==0) ccnc = 2;
00583 if(!(ccnc==1||ccnc==2)) {
00584 MSG("NeugenWC",Msg::kError) << "Invalid event:iaction value "
00585 << iaction << endl;
00586 return 1;
00587 }
00588
00589 int process = 0;
00590 if(ires>=1001&&ires<=1004) process = ires - 1000;
00591 else if(ires>=1&&ires<=4) process = ires;
00592 else {
00593 if(ires==1005) MSG("NeugenWC",Msg::kInfo)
00594 << "Event:iresonance value is 1005 "
00595 << "- no reweighting performed" << endl;
00596 else MSG("NeugenWC",Msg::kError)
00597 << "Invalid event:iresonance (or event:process) value "
00598 << ires << endl;
00599 return 1;
00600 }
00601
00602
00603 if(process==1) {
00604 kid1=1;
00605 kid2=0;
00606 if(kval[1]==999999) {
00607 MSG("NeugenWC",Msg::kWarning)
00608 << "Current event is QEL: q^2 kinematic variable not set!"
00609 << " Returning weight of 1." << endl;
00610 return 1;
00611 }
00612 }
00613 else if(process==2) {
00614 kid1=1;
00615 kid2=2;
00616 if(kval[1]==999999) {
00617 MSG("NeugenWC",Msg::kWarning)
00618 << "Current event is RES: q^2 kinematic variable not set!"
00619 << " Returning weight of 1." << endl;
00620 return 1;
00621 }
00622 if(kval[2]==999999) {
00623 MSG("NeugenWC",Msg::kWarning)
00624 << "Current event is RES: w^2 kinematic variable not set!"
00625 << " Returning weight of 1." << endl;
00626 return 1;
00627 }
00628
00629 if(kval[2]>0) kval[2]=sqrt(kval[2]);
00630 else {
00631 MSG("NeugenWC",Msg::kWarning)
00632 << "w^2 kinematic variable is negative!"
00633 << " Returning weight of 1." << endl;
00634 return 1;
00635 }
00636 }
00637 else if(process==3) {
00638 kid1=3;
00639 kid2=4;
00640 if(kval[3]==999999) {
00641 MSG("NeugenWC",Msg::kWarning)
00642 << "Current event is DIS: x kinematic variable not set!"
00643 << " Returning weight of 1." << endl;
00644 return 1;
00645 }
00646 if(kval[4]==999999) {
00647 MSG("NeugenWC",Msg::kWarning)
00648 << "Current event is DIS: y kinematic variable not set!"
00649 << " Returning weight of 1." << endl;
00650 return 1;
00651 }
00652 }
00653 else if(process==4) {
00654
00655 return 1;
00656 }
00657
00658
00659 TLorentzVector *nu_lv = new TLorentzVector(nu_px,nu_py,nu_pz,nu_en);
00660 TLorentzVector *tar_lv = new TLorentzVector(tar_px,tar_py,tar_pz,tar_en);
00661 interaction *iact = new interaction(flavor_enum(flavor),
00662 nucleus_enum(nucleus),
00663 ccnc_enum(ccnc),
00664 init_state_enum(initial_state));
00665 final_state *fs = new final_state(had_fs);
00666
00667
00668 double weight = 0;
00669 if(event->useStdXSec) {
00670
00671 if(event->stdXSec==0) event->stdXSec =
00672 fNuWrap->offshell_diff_xsec(nu_lv,tar_lv,
00673 iact,process_enum(process),fs,
00674 kinematic_variable_enum(kid1),float(kval[kid1]),
00675 kinematic_variable_enum(kid2),float(kval[kid2]));
00676
00677 if(!fIsOptimized){
00678 weight = fNuWrap->reweight(nu_lv,tar_lv,
00679 kinematic_variable_enum(kid1),float(kval[kid1]),
00680 kinematic_variable_enum(kid2),float(kval[kid2]),
00681 iact,process_enum(process),fs,
00682 fRwtNeuConfig,event->stdXSec, fUpdateConfig);
00683 fUpdateConfig = false;
00684 }else{
00685
00686 weight = fNuWrap->reweight(nu_lv,tar_lv,
00687 kinematic_variable_enum(kid1),float(kval[kid1]),
00688 kinematic_variable_enum(kid2),float(kval[kid2]),
00689 iact,process_enum(process),fs,
00690 fRwtNeuConfig,0,event->stdXSec);
00691 }
00692 }
00693 else {
00694 weight = fNuWrap->reweight(nu_lv,tar_lv,
00695 kinematic_variable_enum(kid1),float(kval[kid1]),
00696 kinematic_variable_enum(kid2),float(kval[kid2]),
00697 iact,process_enum(process),fs,
00698 fRwtNeuConfig);
00699 }
00700
00701 delete nu_lv;
00702 delete tar_lv;
00703 delete iact;
00704 delete fs;
00705 return weight;
00706
00707 }
00708
00709
00710 void NeugenWeightCalculator::PrintReweightConfig(ostream & stream)
00711 {
00712 if(fRwtNeuConfig) fRwtNeuConfig->print(stream);
00713 }
00714 #endif