NuTransSME Class Reference

#include <NuTransSME.h>

List of all members.

Public Member Functions

 NuTransSME (const NuXMLConfig *xmlConfig)
 ~NuTransSME ()
void Fill (NuMCEvent &mc)
void Reweight (NuEvent &nu)
void DoSystematicShifts (NuEvent &event)
void DoSystematicShifts (NuEvent &event, const NuXMLConfig *xmlConfig)
NuEventGetTransitionEvent (const NuEvent &nu)
Bool_t IsSetToDoTransitions () const

Private Member Functions

Float_t Baseline () const
Float_t OscillationWeight (const Float_t energy, const Int_t inu, const NuXMLConfig *lxmlConfig=0) const
Float_t OscillationWeightTrans (const Float_t energy, const Int_t inu, const NuXMLConfig *lxmlConfig=0) const
Double_t GetXSecGraph (NuEvent &nu, bool opp_charge=false) const
Double_t GetXSecGraph (NuMCEvent &nu, bool opp_charge=false) const
Double_t GetXSecNeugen (NuEvent &nu, bool opp_charge=false) const

Private Attributes

TH1D * numu_flux
TH1D * nubar_flux
TH1D * flux_weight_bar
TH1D * flux_weight_nu
TGraph * fNuMuXGraph
TGraph * fNuBarXGraph
NuXMLConfigfxmlConfig
NuSystematicfSyst
bool applyTranSyst
bool firstRunTransBar
bool firstRunTransNu
bool doNothing

Detailed Description

Definition at line 18 of file NuTransSME.h.


Constructor & Destructor Documentation

NuTransSME::NuTransSME ( const NuXMLConfig xmlConfig  ) 

Definition at line 49 of file NuTransSME.cxx.

References NuXMLConfig::AnaVersionString(), applyTranSyst, NuXMLConfig::BinningScheme(), NuXMLConfig::CMuMu(), NuXMLConfig::CTauTau(), NuXMLConfig::DM2Bar(), NuXMLConfig::DM2Nu(), doNothing, firstRunTransBar, firstRunTransNu, flux_weight_bar, flux_weight_nu, NuSystematic::FluxSyst(), fNuBarXGraph, fNuMuXGraph, fxmlConfig, NuXMLConfig::GMu(), NuXMLConfig::GTau(), Msg::kError, Msg::kInfo, MAXMSG, MSG, NuXMLConfig::Name(), nubar_flux, numu_flux, plots(), NuXMLConfig::SN2Bar(), and NuXMLConfig::SN2Nu().

00050 {
00051     // Get our local xmlConfig copy, if it is safe to do so
00052     if (!xmlConfig) {
00053         // If the XML file does not exist, we can't do anything
00054         doNothing = true;
00055         
00056         MAXMSG("NuTransSME",Msg::kInfo,5)
00057         << "No NuXMLConfig object supplied: "
00058         << "performing no systematic shifts, oscillations, or transitions."
00059         << endl;
00060         return;
00061     }
00062     
00063     fxmlConfig = (NuXMLConfig*)xmlConfig->Clone();
00064     firstRunTransBar = true;
00065     firstRunTransNu = true;
00066 
00067     MSG("NuTransSME",Msg::kInfo) << "Applying parameters: " 
00068     << "dm2nu = " << xmlConfig->DM2Nu()   << ", sn2nu = " << xmlConfig->SN2Nu()
00069     << ", dm2bar = " << xmlConfig->DM2Bar() << ", sn2bar = " << xmlConfig->SN2Bar()
00070     << ", gmu = " << xmlConfig->GMu() << ", gtau = " << xmlConfig->GTau() << ", cmumu = " << xmlConfig->CMuMu() << ", ctautau = " << xmlConfig->CTauTau()<<endl;
00071     
00072     doNothing  = false;
00073     
00074     if (doNothing) return;
00075 
00076     const NuPlots* plots=0;    
00077     NuCutter cutter(xmlConfig->AnaVersionString(), plots);
00078 
00079     //Get the binning scheme set up
00080     const NuUtilities cuts;
00081   
00082     //Binning scheme
00083     NuBinningScheme::NuBinningScheme_t binningScheme =
00084         static_cast<NuBinningScheme::NuBinningScheme_t>(xmlConfig->BinningScheme());
00085   
00086     std::vector<Double_t> vReco = cuts.RecoBins(binningScheme);
00087     std::vector<Double_t> vTrue = cuts.TrueBins(binningScheme);
00088     
00089     Int_t numRecoBins = vReco.size() - 1;
00090     Int_t numTrueBins = vTrue.size() - 1;
00091     
00092     Double_t *recoBinsArray = &(vReco[0]);
00093     Double_t *trueBinsArray = &(vTrue[0]);
00094     
00095     // Creat the reweighting histograms
00096     numu_flux   = new TH1D("numu_flux","numu_flux",numTrueBins,trueBinsArray);
00097     nubar_flux  = new TH1D("nubar_flux","nubar_flux",numRecoBins,recoBinsArray);
00098     flux_weight_bar = new TH1D("flux_weight_bar","flux_weight_bar",numRecoBins,recoBinsArray);
00099     flux_weight_nu = new TH1D("flux_weight_nu","flux_weight_nu",numRecoBins,recoBinsArray);
00100     
00101     // Determine how the systematics are applied
00102 
00103     if (NuSystematic::FluxSyst(*xmlConfig)) { 
00104         MSG("NuTransSME", Msg::kInfo) << "Not applying " << xmlConfig->Name() 
00105         << " to transitions b/c it's a flux systematic." << endl;
00106         applyTranSyst = false;
00107     }
00108     else {
00109         applyTranSyst = true;
00110     }
00111     
00112     // Get the XSec graphs from file
00113     TDirectory *start = gDirectory;
00114     
00115     MSG("NuTransSME",Msg::kInfo) << "Get cross-section graphs from file." << endl;
00116 
00117     const char *fxsecfilename = "/minos/app/richa/Rel2.5/NtupleUtils/data/xsec_minos_modbyrs4_v3_5_0_mk.root";
00118    // const char * fxsecfilename = getenv("xf");
00119     
00120     MSG("NuTransSME",Msg::kInfo) << "Opening xsec file: " << fxsecfilename << endl;
00121     TFile *xsecfile = new TFile(fxsecfilename,"READ");
00122   
00123     // Check to see if this failed
00124     if (!xsecfile->IsOpen()) {
00125       MSG("NuTransSME", Msg::kError) << "\n\n"
00126         << "*****************************************************************************\n"
00127         << "Failed to open cross section file.\n"
00128         << "Currently the cross section file is found through the 'xf' environment variable.\n"
00129         << "Perhaps you need to issue a command such as:\n   "
00130         << "   export xf=\"$SRT_PRIVATE_CONTEXT/NtupleUtils/data/xsec_minos_modbyrs4_v3_5_0_mk.root\"\n"
00131         << "*****************************************************************************\n\n"
00132         << endl;
00133       throw runtime_error("Failed to read crosssection file");
00134     }
00135     
00136     // NuMu
00137     TH1F *fhNuMuCCCrossSections = (TH1F*) xsecfile->Get("h_numu_cc_tot");
00138     Float_t *x = new Float_t[fhNuMuCCCrossSections->GetNbinsX()];
00139     Float_t *y = new Float_t[fhNuMuCCCrossSections->GetNbinsX()];
00140     for(int i=0;i<fhNuMuCCCrossSections->GetNbinsX();i++) {
00141         x[i] = fhNuMuCCCrossSections->GetBinCenter(i+1);
00142         y[i] = fhNuMuCCCrossSections->GetBinContent(i+1);
00143     }
00144     fNuMuXGraph = new TGraph(fhNuMuCCCrossSections->GetNbinsX(),x,y);
00145     if (x) {delete[] x; x = 0;}
00146     if (y) {delete[] y; y = 0;}
00147     
00148     // NuMuBar
00149     TH1F *fhNuBarCCCrossSections = (TH1F*) xsecfile->Get("h_numubar_cc_tot");
00150     x = new Float_t[fhNuBarCCCrossSections->GetNbinsX()];
00151     y = new Float_t[fhNuBarCCCrossSections->GetNbinsX()];
00152     for(int i=0;i<fhNuBarCCCrossSections->GetNbinsX();i++) {
00153         x[i] = fhNuBarCCCrossSections->GetBinCenter(i+1);
00154         y[i] = fhNuBarCCCrossSections->GetBinContent(i+1);
00155     }
00156     fNuBarXGraph = new TGraph(fhNuBarCCCrossSections->GetNbinsX(),x,y);
00157     if (x) {delete[] x; x = 0;}
00158     if (y) {delete[] y; y = 0;}
00159     
00160     xsecfile->Close();
00161     if (xsecfile){delete xsecfile; xsecfile = 0;}
00162     
00163     gDirectory = start;
00164 
00165 }

NuTransSME::~NuTransSME (  ) 

Definition at line 168 of file NuTransSME.cxx.

References flux_weight_bar, flux_weight_nu, nubar_flux, and numu_flux.

00169 {
00170     if (numu_flux) delete numu_flux; numu_flux = 0;
00171     if (nubar_flux) delete nubar_flux; nubar_flux = 0;
00172     if (flux_weight_bar) delete flux_weight_bar; flux_weight_bar = 0;
00173     if (flux_weight_nu) delete flux_weight_nu; flux_weight_nu = 0;
00174 }


Member Function Documentation

Float_t NuTransSME::Baseline (  )  const [inline, private]

Definition at line 47 of file NuTransSME.h.

00047 {return 735.0;};

void NuTransSME::DoSystematicShifts ( NuEvent event,
const NuXMLConfig xmlConfig 
)

Definition at line 352 of file NuTransSME.cxx.

References NuXMLConfig::AnaVersionString(), NuXMLConfig::FullTitle(), Msg::kInfo, MAXMSG, NuXMLConfig::Name(), and plots().

00354 {
00355   if (!xmlConfig){
00356     MAXMSG("NuTransSME",Msg::kInfo,5)
00357       << "No NuXMLConfig object supplied: "
00358       << "performing no systematic shifts."
00359       << endl;
00360     return;
00361   } else if (xmlConfig->Name().IsNull()) {
00362     MAXMSG("NuTransSME",Msg::kInfo,5)
00363       << "No Shift supplied in XML. Not doing shifts."
00364       << endl;
00365     return;
00366   } else {
00367     MAXMSG("NuTransSME",Msg::kInfo,5)
00368       << "Performing systematic shift "
00369       << xmlConfig->FullTitle()
00370       << endl;
00371     const NuPlots* plots=0;
00372  
00373     NuCutter cutter(xmlConfig->AnaVersionString(), plots);
00374     static NuSystematic nuSyst(*xmlConfig);
00375     nuSyst.SetNuBarSelector(cutter.GetCut());
00376     nuSyst.SetCCSelector(cutter.GetCut());
00377 
00378     nuSyst.Shift(event);
00379     return;
00380   }
00381 }

void NuTransSME::DoSystematicShifts ( NuEvent event  ) 

Definition at line 384 of file NuTransSME.cxx.

References fSyst, NuXMLConfig::FullTitle(), fxmlConfig, Msg::kInfo, MAXMSG, and NuSystematic::Shift().

Referenced by Reweight().

00385 {
00386     MAXMSG("NuTransSME",Msg::kInfo,5)
00387     << "Performing systematic shift "
00388     << fxmlConfig->FullTitle()
00389     << endl;
00390     fSyst->Shift(event);
00391 }

void NuTransSME::Fill ( NuMCEvent mc  ) 

Definition at line 177 of file NuTransSME.cxx.

References doNothing, GetXSecGraph(), NuMCEvent::iaction, NuMCEvent::inu, Msg::kInfo, MAXMSG, NuMCEvent::neuEnMC, nubar_flux, numu_flux, and NuMCEvent::rw.

Referenced by NuDSTAna::MMTransSME().

00178 {
00179     if (doNothing) return;
00180     
00181     // CC-only
00182     if (mc.iaction==1) { 
00183         double fw = 1;
00184         
00185         // Correct for Cross section
00186         fw *= this->GetXSecGraph(mc);
00187         
00188         if (mc.inu == 14) {
00189             numu_flux->Fill(mc.neuEnMC, mc.rw/fw);
00190             MAXMSG("NuTransSME",Msg::kInfo,5) << "===> Fill numu, xsec=" << fw << endl;
00191         }
00192         else if (mc.inu == -14) {
00193             nubar_flux->Fill(mc.neuEnMC, mc.rw/fw);
00194             MAXMSG("NuTransSME",Msg::kInfo,5) << "===> Fill nubar, xsec=" << fw << endl;
00195         }
00196         else{
00197             MAXMSG("NuTransSME",Msg::kInfo,5) << "===> Fill nothing, ntyp=" << mc.inu << endl;
00198         }
00199     }
00200 
00201 }

NuEvent * NuTransSME::GetTransitionEvent ( const NuEvent nu  ) 

Definition at line 342 of file NuTransSME.cxx.

00343 {
00344 // Adding this because it's declared. It existed only in NuTransitions I think.
00345 // Not sure what it was meant to be. Maybe it can be removed?
00346 // Should ask Richa at some point. - by jcoelho, 09/02/2014
00347              
00348   return 0;
00349 }       

Double_t NuTransSME::GetXSecGraph ( NuMCEvent nu,
bool  opp_charge = false 
) const [private]

Definition at line 510 of file NuTransSME.cxx.

References fNuBarXGraph, fNuMuXGraph, NuMCEvent::inu, NuMCEvent::neuEnMC, and NueConvention::numu.

00511 {
00512     Double_t nuBarXSec = fNuBarXGraph->Eval(nu.neuEnMC,0,"");
00513     Double_t nuMuXSec =  fNuMuXGraph->Eval( nu.neuEnMC,0,"");
00514     
00515     bool numu = (nu.inu > 0);
00516     if (opp_charge) numu = !numu;
00517        
00518     if (numu) return nuMuXSec;
00519     else      return nuBarXSec;
00520 }

Double_t NuTransSME::GetXSecGraph ( NuEvent nu,
bool  opp_charge = false 
) const [private]

Definition at line 497 of file NuTransSME.cxx.

References fNuBarXGraph, fNuMuXGraph, NuEvent::inu, NuEvent::neuEnMC, and NueConvention::numu.

Referenced by Fill().

00498 {
00499     Double_t nuBarXSec = fNuBarXGraph->Eval(nu.neuEnMC,0,"");
00500     Double_t nuMuXSec =  fNuMuXGraph->Eval( nu.neuEnMC,0,"");
00501     
00502     bool numu = (nu.inu > 0);
00503     if (opp_charge) numu = !numu;
00504     
00505     if (numu) return nuMuXSec;
00506     else      return nuBarXSec;
00507 }

Double_t NuTransSME::GetXSecNeugen ( NuEvent nu,
bool  opp_charge = false 
) const [private]

Definition at line 524 of file NuTransSME.cxx.

References e_undefined_init_state, e_vA, e_vbA, e_vbN, e_vbn, e_vbp, e_vN, e_vn, e_vp, NuEvent::hadronicFinalStateMC, NuEvent::iaction, init(), NuEvent::initialStateMC, NuEvent::inu, NuEvent::iresonance, Msg::kError, Msg::kInfo, Msg::kWarning, MSG, NuEvent::neuEnMC, NuEvent::neuPxMC, NuEvent::neuPyMC, NuEvent::neuPzMC, NuEvent::nucleusMC, neugen_wrapper::offshell_diff_xsec(), NuEvent::q2MC, neugen_config::set_best_parameters(), NuEvent::tgtEnMC, NuEvent::tgtPxMC, NuEvent::tgtPyMC, NuEvent::tgtPzMC, NuEvent::w2MC, NuEvent::xMC, and NuEvent::yMC.

00525 {
00527     // Get Cross Sections //
00529     
00530     static neugen_wrapper *fNuWrap = 0;
00531     //static Registry std_registry;   
00532     static neugen_config *std_config = 0;
00533     
00534     if (!fNuWrap) {
00535         //std_registry.Set("neugen:config_name","MODBYRS");
00536         //std_registry.Set("neugen:config_no",4);
00537         
00538         //NuSystematic::SetNeugenDefaults(std_registry);
00539         
00540         std_config = new neugen_config("stdNeuConfig");
00541         std_config->set_best_parameters();
00542         //if(config) SetStandardConfig(config);
00543         
00544         fNuWrap = new neugen_wrapper(std_config);
00545     }
00546     
00547     //event kinematic variables:
00548     double kval[5] = {0,0,0,0,0};
00549     int kid1 = -1;
00550     int kid2 = -1;
00551     kval[1] = nu.q2MC;
00552     kval[2] = nu.w2MC; 
00553     kval[3] = nu.xMC;
00554     kval[4] = nu.yMC;
00555     
00556     //translate from minossoft-speak to neugen-speak:
00557     int flavor = 4;
00558     if(abs(nu.inu)==12) flavor = 1;
00559     else if(abs(nu.inu)==14) flavor = 2; 
00560     else if(abs(nu.inu)==16) flavor = 3;
00561     else if(abs(nu.inu)==1 || abs(nu.inu)==2 || abs(nu.inu)==3) flavor = abs(nu.inu);
00562     else {
00563         MSG("NuTransSME",Msg::kError) << "Invalid event:inu value " << nu.inu << endl;
00564         return 0;
00565     }
00566     
00567     int ccnc = nu.iaction;
00568     if(ccnc==0) ccnc = 2;
00569     if(!(ccnc==1||ccnc==2)) {
00570         MSG("NuTransSME",Msg::kError) << "Invalid event:iaction value " 
00571         << nu.iaction << endl;
00572         return 0;
00573     }
00574     
00575     int process = 0;
00576     if(nu.iresonance>=1001&&nu.iresonance<=1004) process = nu.iresonance - 1000;
00577     else if(nu.iresonance>=1&&nu.iresonance<=4) process = nu.iresonance;
00578     else {
00579         if(nu.iresonance==1005) MSG("NuTransSME",Msg::kInfo) 
00580             << "Event:iresonance value is 1005 "
00581             << "- no reweighting performed" << endl;
00582         else MSG("NeugenWC",Msg::kError) 
00583             << "Invalid event:iresonance (or event:process) value "
00584             << nu.iresonance << endl;
00585         return 0;
00586     }
00587     
00588     //extract necessary kinematic variables:
00589     if(process==1) {   //qel
00590         kid1=1;
00591         kid2=0;
00592         if(kval[1]==999999) {
00593             MSG("NuTransSME",Msg::kWarning) 
00594             << "Current event is QEL: q^2 kinematic variable not set!"
00595             << " Returning weight of 1." << endl;
00596             return 0;
00597         }    
00598     }
00599     else if(process==2) { //res
00600         kid1=1;
00601         kid2=2;
00602         if(kval[1]==999999) {
00603             MSG("NuTransSME",Msg::kWarning) 
00604             << "Current event is RES: q^2 kinematic variable not set!"
00605             << " Returning weight of 1." << endl;
00606             return 0;
00607         }
00608         if(kval[2]==999999) {
00609             MSG("NuTransSME",Msg::kWarning) 
00610             << "Current event is RES: w^2 kinematic variable not set!"
00611             << " Returning weight of 1." << endl;
00612             return 0;
00613         }
00614         //need to pass w rather than w^2 to neugen
00615         if(kval[2]>0) kval[2]=sqrt(kval[2]); 
00616         else {
00617             MSG("NuTransSME",Msg::kWarning) 
00618             << "w^2 kinematic variable is negative!"
00619             << " Returning weight of 1." << endl;
00620             return 0;
00621         }
00622     }
00623     else if(process==3) { //dis
00624         kid1=3;
00625         kid2=4;
00626         if(kval[3]==999999) {
00627             MSG("NuTransSME",Msg::kWarning) 
00628             << "Current event is DIS: x kinematic variable not set!"
00629             << " Returning weight of 1." << endl;
00630             return 0;
00631         }
00632         if(kval[4]==999999) {
00633             MSG("NuTransSME",Msg::kWarning) 
00634             << "Current event is DIS: y kinematic variable not set!"
00635             << " Returning weight of 1." << endl;
00636             return 0;
00637         }
00638     }
00639     else if(process==4) { //coh
00640         //reweighting not possible here yet
00641         return 0;
00642     }
00643     
00644     //start making neugen style objects:
00645     TLorentzVector *nu_lv = new TLorentzVector(nu.neuPxMC,nu.neuPyMC,nu.neuPzMC,nu.neuEnMC);
00646     TLorentzVector *tar_lv = new TLorentzVector(nu.tgtPxMC,nu.tgtPyMC,nu.tgtPzMC,nu.tgtEnMC);
00647     
00648     init_state_t init = init_state_enum(nu.initialStateMC);
00649     
00650     if (opp_charge) {
00651         using namespace init_state;
00652         switch(init) {
00653             case e_vN:   init = e_vbN;  break;
00654             case e_vbN:  init = e_vN;   break;
00655                 
00656             case e_vA:   init = e_vbA;  break;
00657             case e_vbA:  init = e_vA;   break;
00658                 
00659             case e_vp:   init = e_vbp;  break;
00660             case e_vbp:  init = e_vp;   break;
00661                 
00662             case e_vn:   init = e_vbn;  break;
00663             case e_vbn:  init = e_vn;   break;
00664                 
00665             case e_undefined_init_state:
00666             default:            
00667                 MSG("NuTransSME",Msg::kWarning) << "Unknown initial state" << endl;
00668                 return 0;
00669         }
00670         
00671         if (process == 1) {
00672             if (init == e_vbn) init = e_vbp;
00673             if (init == e_vp)  init = e_vn;
00674         }  
00675     }
00676     
00677     interaction *iact = new interaction(flavor_enum(flavor),
00678                                         nucleus_enum(nu.nucleusMC),
00679                                         ccnc_enum(ccnc),
00680                                         init); 
00681     
00682     final_state *fs = new final_state(nu.hadronicFinalStateMC);
00683     
00684     
00685     //calculate cross section
00686     Double_t sigma = fNuWrap->offshell_diff_xsec(nu_lv,tar_lv,
00687                                                  iact,process_enum(process),fs,
00688                                                  kinematic_variable_enum(kid1),float(kval[kid1]),
00689                                                  kinematic_variable_enum(kid2),float(kval[kid2]));
00690     
00691     
00692     delete nu_lv;
00693     delete tar_lv;
00694     delete iact;
00695     delete fs;
00696     
00697     
00698     return sigma;
00699 }

Bool_t NuTransSME::IsSetToDoTransitions (  )  const [inline]

Definition at line 31 of file NuTransSME.h.

00031 {return !doNothing;};

Float_t NuTransSME::OscillationWeight ( const Float_t  energy,
const Int_t  inu,
const NuXMLConfig lxmlConfig = 0 
) const [private]

Definition at line 396 of file NuTransSME.cxx.

References NuXMLConfig::CMuMu(), NuXMLConfig::CTauTau(), NuXMLConfig::DM2Bar(), NuXMLConfig::DM2Nu(), NuXMLConfig::GMu(), NuXMLConfig::GTau(), NuOscProbCalc::OscillationWeightMu2Mu(), NuOscProbCalc::OscillationWeightMu2Tau(), NuXMLConfig::SN2Bar(), and NuXMLConfig::SN2Nu().

Referenced by Reweight().

00397 {
00398     const NuXMLConfig *xmlConfig;
00399     if (lxmlConfig) xmlConfig = lxmlConfig;
00400     else            xmlConfig = fxmlConfig;
00401 
00402     if (14==inu){
00403     return NuOscProbCalc::OscillationWeightMu2Mu(energy,
00404                                           xmlConfig->DM2Nu(),
00405                                           xmlConfig->SN2Nu(),
00406                                           xmlConfig->GMu(),
00407                                           xmlConfig->GTau(),
00408                                           xmlConfig->CMuMu(),
00409                                           xmlConfig->CTauTau());
00410     }   
00411     
00412     else if (-14==inu){                       
00413   
00414     return NuOscProbCalc::OscillationWeightMu2Mu(energy,
00415                                           xmlConfig->DM2Bar(),
00416                                           xmlConfig->SN2Bar(),
00417                                           xmlConfig->GMu(),
00418                                           xmlConfig->GTau(),
00419                                           xmlConfig->CMuMu(),
00420                                           xmlConfig->CTauTau());
00421    }    
00422      
00423    if (16==inu){
00424     return NuOscProbCalc::OscillationWeightMu2Tau(energy,
00425                                             xmlConfig->DM2Nu(),
00426                                             xmlConfig->SN2Nu(),
00427                                             xmlConfig->GMu(),
00428                                             xmlConfig->GTau(),
00429                                             xmlConfig->CMuMu(),
00430                                             xmlConfig->CTauTau());
00431   }
00432   else if (-16==inu){
00433     return NuOscProbCalc::OscillationWeightMu2Tau(energy,
00434                                             xmlConfig->DM2Bar(),
00435                                             xmlConfig->SN2Bar(),
00436                                             xmlConfig->GMu(),
00437                                             xmlConfig->GTau(),
00438                                             xmlConfig->CMuMu(),
00439                                             xmlConfig->CTauTau());
00440   }
00441   
00442     else return 1; // Doesn't oscillate
00443 }

Float_t NuTransSME::OscillationWeightTrans ( const Float_t  energy,
const Int_t  inu,
const NuXMLConfig lxmlConfig = 0 
) const [private]

Definition at line 447 of file NuTransSME.cxx.

References NuXMLConfig::CMuMu(), NuXMLConfig::CTauTau(), NuXMLConfig::DM2Bar(), NuXMLConfig::DM2Nu(), NuXMLConfig::GMu(), NuXMLConfig::GTau(), NuOscProbCalc::OscillationWeightMu2MuBar(), NuOscProbCalc::OscillationWeightMu2TauBar(), NuXMLConfig::SN2Bar(), and NuXMLConfig::SN2Nu().

Referenced by Reweight().

00448 {
00449     const NuXMLConfig *xmlConfig;
00450     if (lxmlConfig) xmlConfig = lxmlConfig;
00451     else            xmlConfig = fxmlConfig;
00452 
00453     if (14==inu){
00454     return NuOscProbCalc::OscillationWeightMu2MuBar(energy,
00455                                           xmlConfig->DM2Nu(),
00456                                           xmlConfig->SN2Nu(),
00457                                           xmlConfig->GMu(),
00458                                           xmlConfig->GTau(),
00459                                           xmlConfig->CMuMu(),
00460                                           xmlConfig->CTauTau());
00461     }
00462     
00463     else if (-14==inu){
00464     
00465     return NuOscProbCalc::OscillationWeightMu2MuBar(energy,
00466                                           xmlConfig->DM2Bar(),
00467                                           xmlConfig->SN2Bar(),
00468                                           xmlConfig->GMu(),
00469                                           xmlConfig->GTau(),
00470                                           xmlConfig->CMuMu(),
00471                                           xmlConfig->CTauTau());
00472    }
00473 
00474    if (16==inu){
00475     return NuOscProbCalc::OscillationWeightMu2TauBar(energy,
00476                                             xmlConfig->DM2Nu(),
00477                                             xmlConfig->SN2Nu(),
00478                                             xmlConfig->GMu(),
00479                                             xmlConfig->GTau(),
00480                                             xmlConfig->CMuMu(),
00481                                             xmlConfig->CTauTau());
00482   }
00483   else if (-16==inu){
00484     return NuOscProbCalc::OscillationWeightMu2TauBar(energy,
00485                                             xmlConfig->DM2Bar(),
00486                                             xmlConfig->SN2Bar(),
00487                                             xmlConfig->GMu(),
00488                                             xmlConfig->GTau(),
00489                                             xmlConfig->CMuMu(),
00490                                             xmlConfig->CTauTau());
00491   }
00492 
00493     else return 1; // Doesn't oscillate
00494 }

void NuTransSME::Reweight ( NuEvent nu  ) 

Definition at line 204 of file NuTransSME.cxx.

References applyTranSyst, NuXMLConfig::CMuMu(), NuXMLConfig::CTauTau(), NuEvent::detector, NuXMLConfig::DM2Bar(), NuXMLConfig::DM2Nu(), doNothing, DoSystematicShifts(), NuEvent::entry, firstRunTransBar, firstRunTransNu, flux_weight_bar, flux_weight_nu, fxmlConfig, NuXMLConfig::GMu(), NuXMLConfig::GTau(), NuEvent::iaction, NuEvent::inu, NuEvent::inunoosc, SimFlag::kData, Detector::kFar, Msg::kInfo, MAXMSG, NuEvent::neuEnMC, nubar_flux, numu_flux, OscillationWeight(), OscillationWeightTrans(), NuEvent::rw, NuEvent::simFlag, NuXMLConfig::SN2Bar(), NuXMLConfig::SN2Nu(), and NuUtilities::ValueAt().

Referenced by NuDSTAna::MMTransSME().

00205 {
00206     // With no xml file, do neither systematics nor oscillations
00207     if (doNothing) return;    
00208     cout<<"____________________________"<<endl;
00209     cout<<"entry="<<nu.entry<<endl;
00210  
00211     double wTotal = 1;
00212 
00213     cout<<"   "<<endl;
00214     // If Systematics are applied to transitioned events,
00215     // can simply shift the whole event
00216 
00217     // nu    -> nubar (unshifted)
00218     // nubar -> nubar (shifted)
00219 
00220     if (applyTranSyst) {
00221         this->DoSystematicShifts(nu,fxmlConfig);
00222         cout<<"doing trans syst"<<endl;
00223     }
00224     
00225     else {
00226         NuEvent * nuCopy = (NuEvent*)nu.Clone();
00227         this->DoSystematicShifts(*nuCopy,fxmlConfig);
00228         wTotal *= nuCopy->rw / nu.rw; 
00229         
00230         cout<<"wTotal #1 ="<<wTotal<<", nuCopy->rw="<<nuCopy->rw<<", nu.rw="<<nu.rw<<endl; 
00231         delete nuCopy;
00232 
00233     }
00234     
00235     //bool isNuE = (nu.inu == 12 || nu.inu == -12);
00236     bool isNuMu = (nu.inu == 14 || nu.inu == -14);
00237     bool isNuTau = (nu.inu == 16 || nu.inu == -16);
00238   
00239   
00240     if (SimFlag::kData == nu.simFlag) {
00241         // don't oscillate real data
00242         MAXMSG("NuTransSME",Msg::kInfo,1) 
00243         << "Not applying fake oscillations to data" << endl;
00244     }
00245     else if (Detector::kFar != nu.detector) { 
00246         // don't oscillate the ND data
00247         MAXMSG("NuTransSME",Msg::kInfo,1) 
00248         << "Not applying fake oscillations to ND MC" << endl;
00249     }
00250     else if (1 != nu.iaction) { 
00251         // check for !CC (i.e. NC)
00252         if(isNuTau) {
00253             MAXMSG("NuTransSME",Msg::kInfo,1) << "Removing NC events from tau file (nu.rw set to zero)" << endl;
00254             wTotal = 0;
00255         }
00256         MAXMSG("NuTransSME",Msg::kInfo,1) << "Not oscillating NCs" << endl;
00257     }
00258     else if ( !(isNuMu || isNuTau) ) { 
00259         // select everything that is not NuMu(bar) or NuTau(bar) i.e. true if nue(bar)
00260         MAXMSG("NuTransSME",Msg::kInfo,1)
00261         << "Not oscillating nue events (non-muon/tau (anti)neutrinos)" << endl;
00262     }
00263     else if ( isNuTau && (nu.inunoosc == 12 || nu.inunoosc==-12) ) {
00264         // select events that started as nues and were turned in to nutaus
00265         MAXMSG("NuTransSME",Msg::kInfo,1)
00266         << "Removing beam nue events from tau file (nu.rw set to zero)" << endl;
00267         wTotal = 0;
00268     }
00269     else if (fxmlConfig->DM2Nu() < 0.0 ||  // Do this last, so that taus and nue's are still properly
00270              fxmlConfig->SN2Nu() < 0.0 ||  // reweighted, even in the no-oscillation case
00271              fxmlConfig->DM2Bar() < 0.0 || // For example, in the no-osc case the nutaus are given a weight of 0
00272              fxmlConfig->SN2Bar() < 0.0){  // so have to explictly just not touch them
00273         MAXMSG("NuTransSME",Msg::kInfo,1)
00274         <<"Not applying fake oscillations due to xml configuration"<<endl;
00275     }
00276         //include the oscillation weight in the general nu.rw variable
00277         
00278     else {
00279         wTotal *= this->OscillationWeight(nu.neuEnMC, nu.inu);
00280         MAXMSG("NuTransSME",Msg::kInfo,10)
00281         << "Oscillating with dm2=" << fxmlConfig->DM2Bar()
00282         << "; weight=" << this->OscillationWeight(nu.neuEnMC, nu.inu)
00283         << "; energy=" << nu.neuEnMC << endl;
00284 
00285         // Only create transitioned events if there are transitions
00286         if (fxmlConfig->GMu() != 0 || fxmlConfig->GTau() != 0 || fxmlConfig->CTauTau() != 0 || fxmlConfig->CMuMu() !=0) {
00287             if (nu.inu == -14) { // and this is an antineutrino
00288                 MAXMSG("NuTransSME",Msg::kInfo,10) << "Getting a transitioned event" << endl;
00289 
00290                 cout<<"Transition Prob = "<<this->OscillationWeightTrans(nu.neuEnMC, nu.inu)<<endl;
00291                 
00292                 if (firstRunTransBar) {
00293                     flux_weight_bar->Divide(numu_flux, nubar_flux);
00294                     firstRunTransBar = false;
00295                 }            
00296                 
00297                 Double_t transitionWeight =  this->OscillationWeightTrans(nu.neuEnMC,nu.inu); 
00298                 //cout<<"Oscillation Weight="<<this->OscillationWeight(nu.neuEnMC,nu.inu)<<endl;
00299                 transitionWeight *= NuUtilities::ValueAt(flux_weight_bar, nu.neuEnMC);
00300                 //cout<<"NuUtilities::ValueAt(flux_weight_bar, nu.neuEnMC)"<<NuUtilities::ValueAt(flux_weight_bar, nu.neuEnMC)<<endl;                
00301                 // Add, not multiply, the transition weight
00302 
00303                 //cout<<"wTotal before adding transwt = "<<wTotal<<endl;
00304                 wTotal += transitionWeight;
00305                 cout<<"inu=="<<nu.inu<<", weight transitioned ="<<wTotal<<endl;
00306  
00307              }
00308 
00309             else if (nu.inu == 14) { // and this is a neutrino
00310                 MAXMSG("NuTransSME",Msg::kInfo,10) << "Getting a transitioned neutrino" << endl;
00311         
00312                 cout<<"Transition Prob Nu = "<<this->OscillationWeightTrans(nu.neuEnMC, nu.inu)<<endl;
00313                 
00314                 if (firstRunTransNu) {
00315                     flux_weight_nu->Divide(nubar_flux, numu_flux);
00316                     firstRunTransNu = false; 
00317                 }            
00318                 
00319                 Double_t transitionWeight =  this->OscillationWeightTrans(nu.neuEnMC,nu.inu);
00320                 //cout<<"Oscillation Weight="<<this->OscillationWeight(nu.neuEnMC,nu.inu)<<endl;
00321                 transitionWeight *= NuUtilities::ValueAt(flux_weight_nu, nu.neuEnMC);
00322                 //cout<<"NuUtilities::ValueAt(flux_weight_nu, nu.neuEnMC)"<<NuUtilities::ValueAt(flux_weight_nu, nu.neuEnMC)<<endl;                
00323                 // Add, not multiply, the transition weight
00324     
00325                 //cout<<"wTotal before adding transwt = "<<wTotal<<endl;
00326                 wTotal += transitionWeight;
00327                 cout<<"weight transitioned for nus="<<wTotal<<endl;
00328         
00329              }
00330 
00331          } //transitions for b or c !=0
00332      } // oscillated everything
00333 
00334     cout<<"wTotal org="<<wTotal<<endl;
00335     nu.rw *= wTotal;
00336     cout<<"inu="<<nu.inu<<", Total nu.rw"<<nu.rw<<endl;
00337   
00338     cout<<"____________________________"<<endl;
00339 }


Member Data Documentation

bool NuTransSME::applyTranSyst [private]

Definition at line 62 of file NuTransSME.h.

Referenced by NuTransSME(), and Reweight().

bool NuTransSME::doNothing [private]

Definition at line 66 of file NuTransSME.h.

Referenced by Fill(), NuTransSME(), and Reweight().

Definition at line 64 of file NuTransSME.h.

Referenced by NuTransSME(), and Reweight().

Definition at line 65 of file NuTransSME.h.

Referenced by NuTransSME(), and Reweight().

TH1D* NuTransSME::flux_weight_bar [private]

Definition at line 37 of file NuTransSME.h.

Referenced by NuTransSME(), Reweight(), and ~NuTransSME().

TH1D* NuTransSME::flux_weight_nu [private]

Definition at line 38 of file NuTransSME.h.

Referenced by NuTransSME(), Reweight(), and ~NuTransSME().

TGraph* NuTransSME::fNuBarXGraph [private]

Definition at line 40 of file NuTransSME.h.

Referenced by GetXSecGraph(), and NuTransSME().

TGraph* NuTransSME::fNuMuXGraph [private]

Definition at line 39 of file NuTransSME.h.

Referenced by GetXSecGraph(), and NuTransSME().

Definition at line 42 of file NuTransSME.h.

Referenced by DoSystematicShifts().

Definition at line 41 of file NuTransSME.h.

Referenced by DoSystematicShifts(), NuTransSME(), and Reweight().

TH1D* NuTransSME::nubar_flux [private]

Definition at line 36 of file NuTransSME.h.

Referenced by Fill(), NuTransSME(), Reweight(), and ~NuTransSME().

TH1D* NuTransSME::numu_flux [private]

Definition at line 31 of file NuTransSME.h.

Referenced by Fill(), NuTransSME(), Reweight(), and ~NuTransSME().


The documentation for this class was generated from the following files:

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1