NuTransition Class Reference

#include <NuTransition.h>

List of all members.

Public Member Functions

 NuTransition (const NuXMLConfig *xmlConfig)
 ~NuTransition ()
void Fill (NuMCEvent &mc)
void Reweight (NuEvent &nu)
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
void DoSystematicShifts (NuEvent &event)
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
TGraph * fNuMuXGraph
TGraph * fNuBarXGraph
NuXMLConfigfxmlConfig
NuSystematicfSyst
bool applyTranSyst
bool firstRunTrans
bool doNothing


Detailed Description

Definition at line 18 of file NuTransition.h.


Constructor & Destructor Documentation

NuTransition::NuTransition ( const NuXMLConfig xmlConfig  ) 

Definition at line 46 of file NuTransition.cxx.

References applyTranSyst, NuXMLConfig::BinningScheme(), NuXMLConfig::DM2Bar(), NuXMLConfig::DM2Nu(), doNothing, firstRunTrans, flux_weight, NuSystematic::FluxSyst(), fNuBarXGraph, fNuMuXGraph, fSyst, fxmlConfig, Msg::kError, Msg::kInfo, MAXMSG, MSG, nubar_flux, numu_flux, NuUtilities::RecoBins(), NuXMLConfig::SN2Bar(), NuXMLConfig::SN2Nu(), NuXMLConfig::TransitionProb(), and NuUtilities::TrueBins().

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

NuTransition::~NuTransition (  ) 

Definition at line 160 of file NuTransition.cxx.

References flux_weight, nubar_flux, and numu_flux.

00161 {
00162     if (numu_flux) delete numu_flux; numu_flux = 0;
00163     if (nubar_flux) delete nubar_flux; nubar_flux = 0;
00164     if (flux_weight) delete flux_weight; flux_weight = 0;
00165 }


Member Function Documentation

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

Definition at line 39 of file NuTransition.h.

00039 {return 735.0;};

void NuTransition::DoSystematicShifts ( NuEvent event  )  [private]

Definition at line 337 of file NuTransition.cxx.

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

Referenced by Reweight().

00338 {
00339     MAXMSG("NuTransition",Msg::kInfo,5)
00340     << "Performing systematic shift "
00341     << fxmlConfig->FullTitle()
00342     << endl;
00343     fSyst->Shift(event);
00344 }

void NuTransition::Fill ( NuMCEvent mc  ) 

Definition at line 168 of file NuTransition.cxx.

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

00169 {
00170     if (doNothing) return;
00171     
00172     // CC-only
00173     if (mc.iaction==1) { 
00174         double fw = 1;
00175         
00176         // Correct for Cross section
00177         fw *= this->GetXSecGraph(mc);
00178         
00179         if (mc.inu == 14) {
00180             numu_flux->Fill(mc.neuEnMC, mc.rw/fw);
00181             MAXMSG("NuTransition",Msg::kInfo,5) << "===> Fill numu, xsec=" << fw << endl;
00182         }
00183         else if (mc.inu == -14) {
00184             nubar_flux->Fill(mc.neuEnMC, mc.rw/fw);
00185             MAXMSG("NuTransition",Msg::kInfo,5) << "===> Fill nubar, xsec=" << fw << endl;
00186         }
00187         else{
00188             MAXMSG("NuTransition",Msg::kInfo,5) << "===> Fill nothing, ntyp=" << mc.inu << endl;
00189         }
00190     }
00191 }

NuEvent * NuTransition::GetTransitionEvent ( const NuEvent nu  ) 

Definition at line 294 of file NuTransition.cxx.

00295 {
00296 // Adding this back with no implementation because it's declared. 
00297 // Not sure what it was meant to be. Maybe it can be removed?
00298 // Should ask Richa at some point. - by jcoelho, 09/02/2014
00299 
00300 /*
00301   MAXMSG("NuTransition",Msg::kInfo,10) << "Getting a transitioned event" << endl;
00302   if (firstRunTrans) {
00303     flux_weight->Divide(numu_flux, nubar_flux);
00304     firstRunTrans = false;
00305   }
00306   
00307   
00308   // Does the event transition?  If not, return an empty pointer.
00309   if (nu.inu != -14 || nu.iaction != 1 ||
00310       fxmlConfig->TransitionProb() <= 0 || nu.detector != Detector::kFar){
00311     return 0;
00312   }
00313   
00314   // Make sure the flux weights have been calculated
00315   
00316   NuEvent *nuTrans = (NuEvent*)nu.Clone();
00317   
00318   if (applyTranSyst) {
00319     this->DoSystematicShifts(*nuTrans);
00320   }
00321   
00322   // Disappearance Probability = 1 - survival probability
00323   Double_t transitionWeight =  1.0 - this->OscillationWeight(nuTrans->neuEnMC,nuTrans->inu,fxmlConfig);
00324   transitionWeight *= fxmlConfig->TransitionProb();
00325   transitionWeight *= NuUtilities::ValueAt(flux_weight, nuTrans->neuEnMC);
00326   
00327   nuTrans->rw *= (transitionWeight);
00328   
00329   if (nuTrans->rw == 0) MSG("NuTransition",Msg::kWarning) << "Warning: zero weight" << endl;
00330   
00331   return nuTrans;
00332 */
00333   return 0;
00334 }           

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

Definition at line 393 of file NuTransition.cxx.

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

00394 {
00395     Double_t nuBarXSec = fNuBarXGraph->Eval(nu.neuEnMC,0,"");
00396     Double_t nuMuXSec =  fNuMuXGraph->Eval( nu.neuEnMC,0,"");
00397     
00398     bool numu = (nu.inu > 0);
00399     if (opp_charge) numu = !numu;
00400        
00401     if (numu) return nuMuXSec;
00402     else      return nuBarXSec;
00403 }

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

Definition at line 380 of file NuTransition.cxx.

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

Referenced by Fill().

00381 {
00382     Double_t nuBarXSec = fNuBarXGraph->Eval(nu.neuEnMC,0,"");
00383     Double_t nuMuXSec =  fNuMuXGraph->Eval( nu.neuEnMC,0,"");
00384     
00385     bool numu = (nu.inu > 0);
00386     if (opp_charge) numu = !numu;
00387     
00388     if (numu) return nuMuXSec;
00389     else      return nuBarXSec;
00390 }

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

Definition at line 407 of file NuTransition.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.

00408 {
00410     // Get Cross Sections //
00412     
00413     static neugen_wrapper *fNuWrap = 0;
00414     //static Registry std_registry;   
00415     static neugen_config *std_config = 0;
00416     
00417     if (!fNuWrap) {
00418         //std_registry.Set("neugen:config_name","MODBYRS");
00419         //std_registry.Set("neugen:config_no",4);
00420         
00421         //NuSystematic::SetNeugenDefaults(std_registry);
00422         
00423         std_config = new neugen_config("stdNeuConfig");
00424         std_config->set_best_parameters();
00425         //if(config) SetStandardConfig(config);
00426         
00427         fNuWrap = new neugen_wrapper(std_config);
00428     }
00429     
00430     //event kinematic variables:
00431     double kval[5] = {0,0,0,0,0};
00432     int kid1 = -1;
00433     int kid2 = -1;
00434     kval[1] = nu.q2MC;
00435     kval[2] = nu.w2MC; 
00436     kval[3] = nu.xMC;
00437     kval[4] = nu.yMC;
00438     
00439     //translate from minossoft-speak to neugen-speak:
00440     int flavor = 4;
00441     if(abs(nu.inu)==12) flavor = 1;
00442     else if(abs(nu.inu)==14) flavor = 2; 
00443     else if(abs(nu.inu)==16) flavor = 3;
00444     else if(abs(nu.inu)==1 || abs(nu.inu)==2 || abs(nu.inu)==3) flavor = abs(nu.inu);
00445     else {
00446         MSG("NuTransition",Msg::kError) << "Invalid event:inu value " << nu.inu << endl;
00447         return 0;
00448     }
00449     
00450     int ccnc = nu.iaction;
00451     if(ccnc==0) ccnc = 2;
00452     if(!(ccnc==1||ccnc==2)) {
00453         MSG("NuTransition",Msg::kError) << "Invalid event:iaction value " 
00454         << nu.iaction << endl;
00455         return 0;
00456     }
00457     
00458     int process = 0;
00459     if(nu.iresonance>=1001&&nu.iresonance<=1004) process = nu.iresonance - 1000;
00460     else if(nu.iresonance>=1&&nu.iresonance<=4) process = nu.iresonance;
00461     else {
00462         if(nu.iresonance==1005) MSG("NuTransition",Msg::kInfo) 
00463             << "Event:iresonance value is 1005 "
00464             << "- no reweighting performed" << endl;
00465         else MSG("NeugenWC",Msg::kError) 
00466             << "Invalid event:iresonance (or event:process) value "
00467             << nu.iresonance << endl;
00468         return 0;
00469     }
00470     
00471     //extract necessary kinematic variables:
00472     if(process==1) {   //qel
00473         kid1=1;
00474         kid2=0;
00475         if(kval[1]==999999) {
00476             MSG("NuTransition",Msg::kWarning) 
00477             << "Current event is QEL: q^2 kinematic variable not set!"
00478             << " Returning weight of 1." << endl;
00479             return 0;
00480         }    
00481     }
00482     else if(process==2) { //res
00483         kid1=1;
00484         kid2=2;
00485         if(kval[1]==999999) {
00486             MSG("NuTransition",Msg::kWarning) 
00487             << "Current event is RES: q^2 kinematic variable not set!"
00488             << " Returning weight of 1." << endl;
00489             return 0;
00490         }
00491         if(kval[2]==999999) {
00492             MSG("NuTransition",Msg::kWarning) 
00493             << "Current event is RES: w^2 kinematic variable not set!"
00494             << " Returning weight of 1." << endl;
00495             return 0;
00496         }
00497         //need to pass w rather than w^2 to neugen
00498         if(kval[2]>0) kval[2]=sqrt(kval[2]); 
00499         else {
00500             MSG("NuTransition",Msg::kWarning) 
00501             << "w^2 kinematic variable is negative!"
00502             << " Returning weight of 1." << endl;
00503             return 0;
00504         }
00505     }
00506     else if(process==3) { //dis
00507         kid1=3;
00508         kid2=4;
00509         if(kval[3]==999999) {
00510             MSG("NuTransition",Msg::kWarning) 
00511             << "Current event is DIS: x kinematic variable not set!"
00512             << " Returning weight of 1." << endl;
00513             return 0;
00514         }
00515         if(kval[4]==999999) {
00516             MSG("NuTransition",Msg::kWarning) 
00517             << "Current event is DIS: y kinematic variable not set!"
00518             << " Returning weight of 1." << endl;
00519             return 0;
00520         }
00521     }
00522     else if(process==4) { //coh
00523         //reweighting not possible here yet
00524         return 0;
00525     }
00526     
00527     //start making neugen style objects:
00528     TLorentzVector *nu_lv = new TLorentzVector(nu.neuPxMC,nu.neuPyMC,nu.neuPzMC,nu.neuEnMC);
00529     TLorentzVector *tar_lv = new TLorentzVector(nu.tgtPxMC,nu.tgtPyMC,nu.tgtPzMC,nu.tgtEnMC);
00530     
00531     init_state_t init = init_state_enum(nu.initialStateMC);
00532     
00533     if (opp_charge) {
00534         using namespace init_state;
00535         switch(init) {
00536             case e_vN:   init = e_vbN;  break;
00537             case e_vbN:  init = e_vN;   break;
00538                 
00539             case e_vA:   init = e_vbA;  break;
00540             case e_vbA:  init = e_vA;   break;
00541                 
00542             case e_vp:   init = e_vbp;  break;
00543             case e_vbp:  init = e_vp;   break;
00544                 
00545             case e_vn:   init = e_vbn;  break;
00546             case e_vbn:  init = e_vn;   break;
00547                 
00548             case e_undefined_init_state:
00549             default:            
00550                 MSG("NuTransition",Msg::kWarning) << "Unknown initial state" << endl;
00551                 return 0;
00552         }
00553         
00554         if (process == 1) {
00555             if (init == e_vbn) init = e_vbp;
00556             if (init == e_vp)  init = e_vn;
00557         }  
00558     }
00559     
00560     interaction *iact = new interaction(flavor_enum(flavor),
00561                                         nucleus_enum(nu.nucleusMC),
00562                                         ccnc_enum(ccnc),
00563                                         init); 
00564     
00565     final_state *fs = new final_state(nu.hadronicFinalStateMC);
00566     
00567     
00568     //calculate cross section
00569     Double_t sigma = fNuWrap->offshell_diff_xsec(nu_lv,tar_lv,
00570                                                  iact,process_enum(process),fs,
00571                                                  kinematic_variable_enum(kid1),float(kval[kid1]),
00572                                                  kinematic_variable_enum(kid2),float(kval[kid2]));
00573     
00574     
00575     delete nu_lv;
00576     delete tar_lv;
00577     delete iact;
00578     delete fs;
00579     
00580     
00581     return sigma;
00582 }

Bool_t NuTransition::IsSetToDoTransitions (  )  const [inline]

Definition at line 27 of file NuTransition.h.

00027 {return !doNothing;};

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

Definition at line 349 of file NuTransition.cxx.

References NuXMLConfig::DM2Bar(), NuXMLConfig::DM2Nu(), NuOscProbCalc::OscillationWeight(), NuXMLConfig::SN2Bar(), and NuXMLConfig::SN2Nu().

Referenced by Reweight().

00350 {
00351     const NuXMLConfig *xmlConfig;
00352     if (lxmlConfig) xmlConfig = lxmlConfig;
00353     else            xmlConfig = fxmlConfig;
00354         
00355     if (14==inu){
00356       return NuOscProbCalc::OscillationWeight(energy,
00357                                             xmlConfig->DM2Nu(),
00358                                             xmlConfig->SN2Nu());
00359     }
00360     else if (-14==inu){
00361       return NuOscProbCalc::OscillationWeight(energy,
00362                                             xmlConfig->DM2Bar(),
00363                                             xmlConfig->SN2Bar());
00364     }
00365     if (16==inu){
00366       return 1-NuOscProbCalc::OscillationWeight(energy,
00367                                               xmlConfig->DM2Nu(),
00368                                               xmlConfig->SN2Nu());
00369     }
00370     else if (-16==inu){
00371       return 1-NuOscProbCalc::OscillationWeight(energy,
00372                                               xmlConfig->DM2Bar(),
00373                                               xmlConfig->SN2Bar());
00374     }
00375     else return 1; // Doesn't oscillate
00376 }

void NuTransition::Reweight ( NuEvent nu  ) 

Definition at line 194 of file NuTransition.cxx.

References applyTranSyst, NuEvent::detector, NuXMLConfig::DM2Bar(), NuXMLConfig::DM2Nu(), doNothing, DoSystematicShifts(), firstRunTrans, flux_weight, fxmlConfig, NuEvent::iaction, NuEvent::inu, NuEvent::inunoosc, SimFlag::kData, Detector::kFar, Msg::kInfo, MAXMSG, NuEvent::neuEnMC, nubar_flux, numu_flux, OscillationWeight(), NuEvent::rw, NuEvent::simFlag, NuXMLConfig::SN2Bar(), NuXMLConfig::SN2Nu(), NuXMLConfig::TransitionProb(), and NuUtilities::ValueAt().

Referenced by DataMCPlots::PrepareEventForPlotting(), and NuDSTAna::StdNMBAna().

00195 {
00196     // With no xml file, do neither systematics or oscillations
00197     if (doNothing) return;    
00198  
00199     double wTotal = 1;
00200     
00201     // If Systematics are applied to transitioned events,
00202     // can simply shift the whole event
00203     if (applyTranSyst) {
00204         this->DoSystematicShifts(nu);
00205     }
00206     else {
00207         NuEvent * nuCopy = (NuEvent*)nu.Clone();
00208         this->DoSystematicShifts(*nuCopy);
00209         wTotal *= nuCopy->rw / nu.rw;
00210         delete nuCopy;
00211     }
00212     
00213   //bool isNuE = (nu.inu == 12 || nu.inu == -12);
00214     bool isNuMu = (nu.inu == 14 || nu.inu == -14);
00215     bool isNuTau = (nu.inu == 16 || nu.inu == -16);
00216   
00217   
00218     if (SimFlag::kData == nu.simFlag) {
00219         // don't oscillate real data
00220         MAXMSG("NuTransition",Msg::kInfo,1) 
00221         << "Not applying fake oscillations to data" << endl;
00222     }
00223     else if (Detector::kFar != nu.detector) { 
00224         // don't oscillate the ND data
00225         MAXMSG("NuTransition",Msg::kInfo,1) 
00226         << "Not applying fake oscillations to ND MC" << endl;
00227     }
00228     else if (1 != nu.iaction) { 
00229         // check for !CC (i.e. NC)
00230         if(isNuTau) {
00231             MAXMSG("NuTransition",Msg::kInfo,1) << "Removing NC events from tau file (nu.rw set to zero)" << endl;
00232             wTotal = 0;
00233         }
00234         MAXMSG("NuTransition",Msg::kInfo,1) << "Not oscillating NCs" << endl;
00235     }
00236     else if ( !(isNuMu || isNuTau) ) { 
00237         // select everything that is not NuMu(bar) or NuTau(bar) i.e. true if nue(bar)
00238         MAXMSG("NuTransition",Msg::kInfo,1)
00239         << "Not oscillating nue events (non-muon/tau (anti)neutinos)" << endl;
00240     }
00241     else if ( isNuTau && (nu.inunoosc == 12 || nu.inunoosc==-12) ) {
00242         // select events that started as nues and were turned in to nutaus
00243         MAXMSG("NuTransition",Msg::kInfo,1)
00244         << "Removing beam nue events from tau file (nu.rw set to zero)" << endl;
00245         wTotal = 0;
00246     }
00247     else if (fxmlConfig->DM2Nu() < 0.0 ||  // Do this last, so that taus and nue's are still properly
00248              fxmlConfig->SN2Nu() < 0.0 ||  // reweighted, even in the no-oscillation case
00249              fxmlConfig->DM2Bar() < 0.0 || // For example, in the no-osc case the nutaus are given a weight of 0
00250              fxmlConfig->SN2Bar() < 0.0){  // so have to explictly just not touch them
00251         MAXMSG("NuTransition",Msg::kInfo,1)
00252         <<"Not applying fake oscillations due to xml configuration"<<endl;
00253     }
00254     else {
00255         //include the oscillation weight in the general nu.rw variable
00256         wTotal *= this->OscillationWeight(nu.neuEnMC, nu.inu);
00257         MAXMSG("NuTransition",Msg::kInfo,10)
00258         << "Oscillating with dm2=" << fxmlConfig->DM2Nu()
00259         << "; weight=" << this->OscillationWeight(nu.neuEnMC, nu.inu)
00260         << "; energy=" << nu.neuEnMC << endl;
00261         
00262         // Only create transitioned events if there are transitions
00263         if (fxmlConfig->TransitionProb() > 0) {
00264             if (isNuTau) {
00265                 MAXMSG("NuTransition",Msg::kInfo,1)
00266                 << "Weighting down appeared NuTaus for transitions by " 
00267                 << (1.0 - fxmlConfig->TransitionProb()) << endl;
00268                 wTotal *= (1.0 - fxmlConfig->TransitionProb());
00269             }
00270             else if (nu.inu == -14) { // and this is an antineutrino
00271                 MAXMSG("NuTransition",Msg::kInfo,10) << "Getting a transitioned event" << endl;
00272                 
00273                 if (firstRunTrans) {
00274                     flux_weight->Divide(numu_flux, nubar_flux);
00275                     firstRunTrans = false;
00276                 }            
00277                 
00278                 // Disappearance Probability = 1 - survival probability
00279                 Double_t transitionWeight =  1.0 - this->OscillationWeight(nu.neuEnMC,nu.inu); 
00280                 transitionWeight *= fxmlConfig->TransitionProb();
00281                 transitionWeight *= NuUtilities::ValueAt(flux_weight, nu.neuEnMC);
00282                 
00283                 // Add, not multiply, the transition weight
00284                 wTotal += transitionWeight;
00285             }
00286         }
00287     }
00288     
00289     nu.rw *= wTotal;
00290 }


Member Data Documentation

bool NuTransition::applyTranSyst [private]

Definition at line 49 of file NuTransition.h.

Referenced by NuTransition(), and Reweight().

bool NuTransition::doNothing [private]

Definition at line 51 of file NuTransition.h.

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

bool NuTransition::firstRunTrans [private]

Definition at line 50 of file NuTransition.h.

Referenced by NuTransition(), and Reweight().

TH1D* NuTransition::flux_weight [private]

Definition at line 33 of file NuTransition.h.

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

TGraph* NuTransition::fNuBarXGraph [private]

Definition at line 35 of file NuTransition.h.

Referenced by GetXSecGraph(), and NuTransition().

TGraph* NuTransition::fNuMuXGraph [private]

Definition at line 34 of file NuTransition.h.

Referenced by GetXSecGraph(), and NuTransition().

NuSystematic* NuTransition::fSyst [private]

Definition at line 37 of file NuTransition.h.

Referenced by DoSystematicShifts(), and NuTransition().

NuXMLConfig* NuTransition::fxmlConfig [private]

Definition at line 36 of file NuTransition.h.

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

TH1D* NuTransition::nubar_flux [private]

Definition at line 32 of file NuTransition.h.

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

TH1D* NuTransition::numu_flux [private]

Definition at line 27 of file NuTransition.h.

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


The documentation for this class was generated from the following files:
Generated on Mon Nov 10 00:56:10 2014 for loon by  doxygen 1.4.7