loon: BMPTWeightCalculator.cxx Source File

00001 #ifndef bmptweightcalculator_cxx
00002 #define bmptweightcalculator_cxx
00003 #include "MessageService/MsgService.h"
00004 #include "MCReweight/BMPTWeightCalculator.h"
00005 #include "MCReweight/BMPTEvent.h"
00006 
00007 CVSID("$Id: BMPTWeightCalculator.cxx,v 1.8 2005/06/02 20:22:19 cbs Exp $");
00008 
00009 //*********************************************
00010 BMPTWeightCalculator::BMPTWeightCalculator(Registry *config) 
00011   : WeightCalculator(config)
00012 {
00013 
00014   fWCname = "BMPTWeightCalculator";
00015 
00016   fStandardConfigChanged = false;
00017   fReweightConfigChanged = false;
00018 
00019   fStdConfig = new BMPTConfig();
00020   if(config) SetStandardConfig(config);
00021   fRwtConfig = new BMPTConfig();
00022   if(config) SetReweightConfig(config);
00023 
00024   fCalc = new BMPTCalc(fStdConfig);
00025   MSG("BMPT",Msg::kDebug) << "In BMPTWeightCalculator Constructor" << endl;  
00026 }
00027 
00028 //*********************************************
00029 BMPTWeightCalculator::~BMPTWeightCalculator()
00030 {
00031   delete fCalc;
00032   delete fStdConfig;
00033   delete fRwtConfig;
00034   MSG("BMPT",Msg::kDebug) << "In BMPTWeightCalculator Destructor" << endl;
00035 }
00036 
00037 //*********************************************
00038 void BMPTWeightCalculator::ReweightConfigReset()
00039 {
00040   fRwtConfig->SetDefaults();
00041   if(fStandardConfig) SetReweightConfig(fStandardConfig);
00042 }
00043 
00044 //*********************************************
00045 void BMPTWeightCalculator::Config()
00046 {
00047 
00048   double tempd = 0;
00049   if(fStandardConfigChanged) {
00050     if(fStandardConfig->Get("bmpt:cpa_pi",tempd)) 
00051       fStdConfig->SetCPA_PI(tempd);
00052     if(fStandardConfig->Get("bmpt:cpa_k",tempd)) 
00053       fStdConfig->SetCPA_K(tempd);
00054     if(fStandardConfig->Get("bmpt:cpa_p",tempd)) 
00055       fStdConfig->SetCPA_P(tempd);
00056     if(fStandardConfig->Get("bmpt:cpa_pbar",tempd)) 
00057       fStdConfig->SetCPA_PBAR(tempd);
00058     if(fStandardConfig->Get("bmpt:cpb_pi",tempd)) 
00059       fStdConfig->SetCPB_PI(tempd);
00060     if(fStandardConfig->Get("bmpt:cpb_k",tempd)) 
00061       fStdConfig->SetCPB_K(tempd);
00062     if(fStandardConfig->Get("bmpt:cpb_p",tempd)) 
00063       fStdConfig->SetCPB_P(tempd);
00064     if(fStandardConfig->Get("bmpt:cpb_pbar",tempd)) 
00065       fStdConfig->SetCPB_PBAR(tempd);
00066     if(fStandardConfig->Get("bmpt:alpha_pi",tempd)) 
00067       fStdConfig->SetALPHA_PI(tempd);
00068     if(fStandardConfig->Get("bmpt:alpha_k",tempd)) 
00069       fStdConfig->SetALPHA_K(tempd);
00070     if(fStandardConfig->Get("bmpt:alpha_p",tempd)) 
00071       fStdConfig->SetALPHA_P(tempd);
00072     if(fStandardConfig->Get("bmpt:alpha_pbar",tempd)) 
00073       fStdConfig->SetALPHA_PBAR(tempd);
00074     if(fStandardConfig->Get("bmpt:beta_pi",tempd)) 
00075       fStdConfig->SetBETA_PI(tempd);
00076     if(fStandardConfig->Get("bmpt:beta_k",tempd)) 
00077       fStdConfig->SetBETA_K(tempd);
00078     if(fStandardConfig->Get("bmpt:beta_p",tempd)) 
00079       fStdConfig->SetBETA_P(tempd);
00080     if(fStandardConfig->Get("bmpt:beta_pbar",tempd)) 
00081       fStdConfig->SetBETA_PBAR(tempd);
00082     if(fStandardConfig->Get("bmpt:gamma_pi",tempd)) 
00083       fStdConfig->SetGAMMA_PI(tempd);
00084     if(fStandardConfig->Get("bmpt:gamma_k",tempd)) 
00085       fStdConfig->SetGAMMA_K(tempd);
00086     if(fStandardConfig->Get("bmpt:gamma_p",tempd)) 
00087       fStdConfig->SetGAMMA_P(tempd);
00088     if(fStandardConfig->Get("bmpt:gamma_pbar",tempd)) 
00089       fStdConfig->SetGAMMA_PBAR(tempd);
00090     if(fStandardConfig->Get("bmpt:delta_pi",tempd)) 
00091       fStdConfig->SetDELTA_PI(tempd);
00092     if(fStandardConfig->Get("bmpt:delta_k",tempd)) 
00093       fStdConfig->SetDELTA_K(tempd);
00094     if(fStandardConfig->Get("bmpt:delta_p",tempd)) 
00095       fStdConfig->SetDELTA_P(tempd);
00096     if(fStandardConfig->Get("bmpt:delta_pbar",tempd)) 
00097       fStdConfig->SetDELTA_PBAR(tempd);
00098     if(fStandardConfig->Get("bmpt:sma_pi",tempd)) 
00099       fStdConfig->SetSMA_PI(tempd);
00100     if(fStandardConfig->Get("bmpt:sma_k",tempd)) 
00101       fStdConfig->SetSMA_K(tempd);
00102     if(fStandardConfig->Get("bmpt:sma_p",tempd)) 
00103       fStdConfig->SetSMA_P(tempd);
00104     if(fStandardConfig->Get("bmpt:sma_pbar",tempd)) 
00105       fStdConfig->SetSMA_PBAR(tempd);
00106     if(fStandardConfig->Get("bmpt:smb_pi",tempd)) 
00107       fStdConfig->SetSMB_PI(tempd);
00108     if(fStandardConfig->Get("bmpt:smb_k",tempd)) 
00109       fStdConfig->SetSMB_K(tempd);
00110     if(fStandardConfig->Get("bmpt:smb_p",tempd)) 
00111       fStdConfig->SetSMB_P(tempd);
00112     if(fStandardConfig->Get("bmpt:smb_pbar",tempd)) 
00113       fStdConfig->SetSMB_PBAR(tempd);
00114     if(fStandardConfig->Get("bmpt:r0_pi",tempd)) 
00115       fStdConfig->SetR0_PI(tempd);
00116     if(fStandardConfig->Get("bmpt:r0_k",tempd)) 
00117       fStdConfig->SetR0_K(tempd);
00118     if(fStandardConfig->Get("bmpt:r0_p",tempd)) 
00119       fStdConfig->SetR0_P(tempd);
00120     if(fStandardConfig->Get("bmpt:r0_pbar",tempd)) 
00121       fStdConfig->SetR0_PBAR(tempd);
00122     if(fStandardConfig->Get("bmpt:r1_pi",tempd)) 
00123       fStdConfig->SetR1_PI(tempd);
00124     if(fStandardConfig->Get("bmpt:r1_k",tempd)) 
00125       fStdConfig->SetR1_K(tempd);
00126     if(fStandardConfig->Get("bmpt:r1_p",tempd)) 
00127       fStdConfig->SetR1_P(tempd);
00128     if(fStandardConfig->Get("bmpt:r1_pbar",tempd)) 
00129       fStdConfig->SetR1_PBAR(tempd);
00130     if(fStandardConfig->Get("bmpt:beam_p",tempd)) 
00131       fStdConfig->SetBEAM_P(tempd);
00132     if(fStandardConfig->Get("bmpt:a_target",tempd)) 
00133       fStdConfig->SetA_TARGET(tempd);
00134     if(fStandardConfig->Get("bmpt:a_be",tempd)) 
00135       fStdConfig->SetA_BE(tempd);
00136     if(fStandardConfig->Get("bmpt:target_rho",tempd)) 
00137       fStdConfig->SetTARGET_RHO(tempd);
00138     if(fStandardConfig->Get("bmpt:lambda_p",tempd)) 
00139       fStdConfig->SetLAMBDA_P(tempd);
00140     if(fStandardConfig->Get("bmpt:lambda_s",tempd)) 
00141       fStdConfig->SetLAMBDA_S(tempd);    
00142     fStandardConfigChanged = false;
00143   }
00144   if(fReweightConfigChanged) {    
00145     //toggle using input values as absolute or scale factors:
00146     Float_t sc = 0; Int_t tempi = 0;
00147     if(fReweightConfig->Get("bmpt:use_scale_factors",tempi)) {
00148       if(tempi==1) sc = 1.0;
00149     }
00150 
00151     //set values:
00152     if(fReweightConfig->Get("bmpt:cpa_pi",tempd))
00153       fRwtConfig->SetCPA_PI(tempd*(1+sc*(fStdConfig->GetCPA_PI()-1))); 
00154     if(fReweightConfig->Get("bmpt:cpa_k",tempd)) 
00155       fRwtConfig->SetCPA_K(tempd*(1+sc*(fStdConfig->GetCPA_K()-1))); 
00156     if(fReweightConfig->Get("bmpt:cpa_p",tempd)) 
00157       fRwtConfig->SetCPA_P(tempd*(1+sc*(fStdConfig->GetCPA_P()-1))); 
00158     if(fReweightConfig->Get("bmpt:cpa_pbar",tempd)) 
00159       fRwtConfig->SetCPA_PBAR(tempd*(1+sc*(fStdConfig->GetCPA_PBAR()-1))); 
00160     if(fReweightConfig->Get("bmpt:cpb_pi",tempd)) 
00161       fRwtConfig->SetCPB_PI(tempd*(1+sc*(fStdConfig->GetCPB_PI()-1))); 
00162     if(fReweightConfig->Get("bmpt:cpb_k",tempd)) 
00163       fRwtConfig->SetCPB_K(tempd*(1+sc*(fStdConfig->GetCPB_K()-1))); 
00164     if(fReweightConfig->Get("bmpt:cpb_p",tempd)) 
00165       fRwtConfig->SetCPB_P(tempd*(1+sc*(fStdConfig->GetCPB_P()-1))); 
00166     if(fReweightConfig->Get("bmpt:cpb_pbar",tempd)) 
00167       fRwtConfig->SetCPB_PBAR(tempd*(1+sc*(fStdConfig->GetCPB_PBAR()-1))); 
00168     if(fReweightConfig->Get("bmpt:alpha_pi",tempd)) 
00169       fRwtConfig->SetALPHA_PI(tempd*(1+sc*(fStdConfig->GetALPHA_PI()-1))); 
00170     if(fReweightConfig->Get("bmpt:alpha_k",tempd)) 
00171       fRwtConfig->SetALPHA_K(tempd*(1+sc*(fStdConfig->GetALPHA_K()-1))); 
00172     if(fReweightConfig->Get("bmpt:alpha_p",tempd)) 
00173       fRwtConfig->SetALPHA_P(tempd*(1+sc*(fStdConfig->GetALPHA_P()-1))); 
00174     if(fReweightConfig->Get("bmpt:alpha_pbar",tempd)) 
00175       fRwtConfig->SetALPHA_PBAR(tempd*(1+sc*(fStdConfig->GetALPHA_PBAR()-1))); 
00176     if(fReweightConfig->Get("bmpt:beta_pi",tempd)) 
00177       fRwtConfig->SetBETA_PI(tempd*(1+sc*(fStdConfig->GetBETA_PI()-1))); 
00178     if(fReweightConfig->Get("bmpt:beta_k",tempd)) 
00179       fRwtConfig->SetBETA_K(tempd*(1+sc*(fStdConfig->GetBETA_K()-1))); 
00180     if(fReweightConfig->Get("bmpt:beta_p",tempd)) 
00181       fRwtConfig->SetBETA_P(tempd*(1+sc*(fStdConfig->GetBETA_P()-1))); 
00182     if(fReweightConfig->Get("bmpt:beta_pbar",tempd)) 
00183       fRwtConfig->SetBETA_PBAR(tempd*(1+sc*(fStdConfig->GetBETA_PBAR()-1))); 
00184     if(fReweightConfig->Get("bmpt:gamma_pi",tempd)) 
00185       fRwtConfig->SetGAMMA_PI(tempd*(1+sc*(fStdConfig->GetGAMMA_PI()-1))); 
00186     if(fReweightConfig->Get("bmpt:gamma_k",tempd)) 
00187       fRwtConfig->SetGAMMA_K(tempd*(1+sc*(fStdConfig->GetGAMMA_K()-1))); 
00188     if(fReweightConfig->Get("bmpt:gamma_p",tempd)) 
00189       fRwtConfig->SetGAMMA_P(tempd*(1+sc*(fStdConfig->GetGAMMA_P()-1))); 
00190     if(fReweightConfig->Get("bmpt:gamma_pbar",tempd)) 
00191       fRwtConfig->SetGAMMA_PBAR(tempd*(1+sc*(fStdConfig->GetGAMMA_PBAR()-1))); 
00192     if(fReweightConfig->Get("bmpt:delta_pi",tempd)) 
00193       fRwtConfig->SetDELTA_PI(tempd*(1+sc*(fStdConfig->GetDELTA_PI()-1))); 
00194     if(fReweightConfig->Get("bmpt:delta_k",tempd)) 
00195       fRwtConfig->SetDELTA_K(tempd*(1+sc*(fStdConfig->GetDELTA_K()-1))); 
00196     if(fReweightConfig->Get("bmpt:delta_p",tempd));
00197       fRwtConfig->SetDELTA_P(tempd*(1+sc*(fStdConfig->GetDELTA_P()-1))); 
00198     if(fReweightConfig->Get("bmpt:delta_pbar",tempd)) 
00199       fRwtConfig->SetDELTA_PBAR(tempd*(1+sc*(fStdConfig->GetDELTA_PBAR()-1))); 
00200     if(fReweightConfig->Get("bmpt:sma_pi",tempd)) 
00201       fRwtConfig->SetSMA_PI(tempd*(1+sc*(fStdConfig->GetSMA_PI()-1))); 
00202     if(fReweightConfig->Get("bmpt:sma_k",tempd)) 
00203       fRwtConfig->SetSMA_K(tempd*(1+sc*(fStdConfig->GetSMA_K()-1))); 
00204     if(fReweightConfig->Get("bmpt:sma_p",tempd)) 
00205       fRwtConfig->SetSMA_P(tempd*(1+sc*(fStdConfig->GetSMA_P()-1))); 
00206     if(fReweightConfig->Get("bmpt:sma_pbar",tempd)) 
00207       fRwtConfig->SetSMA_PBAR(tempd*(1+sc*(fStdConfig->GetSMA_PBAR()-1))); 
00208     if(fReweightConfig->Get("bmpt:smb_pi",tempd)) 
00209       fRwtConfig->SetSMB_PI(tempd*(1+sc*(fStdConfig->GetSMB_PI()-1))); 
00210     if(fReweightConfig->Get("bmpt:smb_k",tempd)) 
00211       fRwtConfig->SetSMB_K(tempd*(1+sc*(fStdConfig->GetSMB_K()-1))); 
00212     if(fReweightConfig->Get("bmpt:smb_p",tempd)) 
00213       fRwtConfig->SetSMB_P(tempd*(1+sc*(fStdConfig->GetSMB_P()-1))); 
00214     if(fReweightConfig->Get("bmpt:smb_pbar",tempd)) 
00215       fRwtConfig->SetSMB_PBAR(tempd*(1+sc*(fStdConfig->GetSMB_PBAR()-1))); 
00216     if(fReweightConfig->Get("bmpt:r0_pi",tempd)) 
00217       fRwtConfig->SetR0_PI(tempd*(1+sc*(fStdConfig->GetR0_PI()-1))); 
00218     if(fReweightConfig->Get("bmpt:r0_k",tempd)) 
00219       fRwtConfig->SetR0_K(tempd*(1+sc*(fStdConfig->GetR0_K()-1))); 
00220     if(fReweightConfig->Get("bmpt:r0_p",tempd)) 
00221       fRwtConfig->SetR0_P(tempd*(1+sc*(fStdConfig->GetR0_P()-1))); 
00222     if(fReweightConfig->Get("bmpt:r0_pbar",tempd)) 
00223       fRwtConfig->SetR0_PBAR(tempd*(1+sc*(fStdConfig->GetR0_PBAR()-1))); 
00224     if(fReweightConfig->Get("bmpt:r1_pi",tempd)) 
00225       fRwtConfig->SetR1_PI(tempd*(1+sc*(fStdConfig->GetR1_PI()-1))); 
00226     if(fReweightConfig->Get("bmpt:r1_k",tempd)) 
00227       fRwtConfig->SetR1_K(tempd*(1+sc*(fStdConfig->GetR1_K()-1))); 
00228     if(fReweightConfig->Get("bmpt:r1_p",tempd)) 
00229       fRwtConfig->SetR1_P(tempd*(1+sc*(fStdConfig->GetR1_P()-1))); 
00230     if(fReweightConfig->Get("bmpt:r1_pbar",tempd)) 
00231       fRwtConfig->SetR1_PBAR(tempd*(1+sc*(fStdConfig->GetR1_PBAR()-1))); 
00232 
00233     //probably shouldn't be trying to change these:
00234     if(fReweightConfig->Get("bmpt:beam_p",tempd)) {
00235       fRwtConfig->SetBEAM_P(tempd*(1+sc*(fStdConfig->GetBEAM_P()-1))); 
00236       MSG("BMPT",Msg::kWarning) << "Changing BEAM_P: are you sure you want to do this?" << endl;  
00237     }
00238     if(fReweightConfig->Get("bmpt:a_target",tempd)) {
00239       fRwtConfig->SetA_TARGET(tempd*(1+sc*(fStdConfig->GetA_TARGET()-1))); 
00240       MSG("BMPT",Msg::kWarning) << "Changing A_TARGET: are you sure you want to do this?" << endl;  
00241     }
00242     if(fReweightConfig->Get("bmpt:a_be",tempd)) {
00243       fRwtConfig->SetA_BE(tempd*(1+sc*(fStdConfig->GetA_BE()-1))); 
00244       MSG("BMPT",Msg::kWarning) << "Changing A_BE: are you sure you want to do this?" << endl;  
00245     }
00246     if(fReweightConfig->Get("bmpt:target_rho",tempd)) {
00247       fRwtConfig->SetTARGET_RHO(tempd*(1+sc*(fStdConfig->GetTARGET_RHO()-1))); 
00248       MSG("BMPT",Msg::kWarning) << "Changing TARGET_RHO: are you sure you want to do this?" << endl;  
00249     }
00250     if(fReweightConfig->Get("bmpt:lambda_p",tempd)) {
00251       fRwtConfig->SetLAMBDA_P(tempd*(1+sc*(fStdConfig->GetLAMBDA_P()-1))); 
00252       MSG("BMPT",Msg::kWarning) << "Changing LAMBDA_P: are you sure you want to do this?" << endl;  
00253     }
00254     if(fReweightConfig->Get("bmpt:lambda_s",tempd)) {
00255       fRwtConfig->SetLAMBDA_S(tempd*(1+sc*(fStdConfig->GetLAMBDA_S()-1))); 
00256       MSG("BMPT",Msg::kWarning) << "Changing LAMBDA_S: are you sure you want to do this?" << endl;  
00257     }
00258     fReweightConfigChanged = false;
00259   }
00260 }
00261 
00262 //*********************************************
00263 double BMPTWeightCalculator::GetWeight(Registry *event)
00264 {
00265 
00266   if(!event) return 1;
00267   //neutrino parent characteristics:
00268   double x = 0;
00269   double y = 0;
00270   double z = 0;
00271   double px = 0;
00272   double py = 0;
00273   double pz = 0;
00274   int pid = 0;
00275   int gen = 0;
00276   if(!event->Get("event:nuparent_x",x)) return 1;
00277   if(!event->Get("event:nuparent_y",y)) return 1;
00278   if(!event->Get("event:nuparent_z",z)) return 1;
00279   if(!event->Get("event:nuparent_px",px)) return 1;
00280   if(!event->Get("event:nuparent_py",py)) return 1;
00281   if(!event->Get("event:nuparent_pz",pz)) return 1;
00282   if(!event->Get("event:nuparent_pid",pid)) return 1;
00283   if(!event->Get("event:nuparent_gen",gen)) return 1;
00284   
00285   BMPTEvent ev(x,y,z,px,py,pz,pid,gen);
00286   return fCalc->Reweight(&ev,fRwtConfig);
00287   
00288 }
00289 
00290 //*********************************************
00291 double BMPTWeightCalculator::GetWeight(MCEventInfo *event,NuParent *parent)
00292 {
00293   
00294   if(!event) {    //surprising but not fatal
00295     MSG("BMPT",Msg::kDebug) << "No MCEventInfo object" << endl;
00296   }
00297 
00298   if(!parent) return 1;
00299   //neutrino parent characteristics:
00300   double x = parent->GetX();
00301   double y = parent->GetY();
00302   double z = parent->GetZ();
00303   double px = parent->GetPx();
00304   double py = parent->GetPy();
00305   double pz = parent->GetPz();
00306   int pid = parent->GetPID();
00307   int gen = parent->GetGen();
00308   if(pz==0) return 1;
00309   
00310   BMPTEvent ev(x,y,z,px,py,pz,pid,gen);
00311   return fCalc->Reweight(&ev,fRwtConfig);
00312   
00313 }
00314 
00315 //*********************************************
00316 void BMPTWeightCalculator::PrintReweightConfig(ostream & stream)
00317 {
00318   if(fRwtConfig) fRwtConfig->Print(stream);
00319 }
00320 #endif