NuTransition.cxx

Go to the documentation of this file.

#include "MessageService/MsgFormat.h"
#include "MessageService/MsgService.h"
#include "Conventions/ReleaseType.h"
#include "Conventions/Munits.h"
#include "Conventions/Detector.h"
#include "Conventions/SimFlag.h"

#include "NuMuBar/NuTransition.h"

#include "NtupleUtils/NuUtilities.h"
#include "NtupleUtils/NuEvent.h"
#include "NtupleUtils/NuMCEvent.h"
#include "NtupleUtils/NuSystematic.h"
#include "NtupleUtils/NuXMLConfig.h"

// Neugen Stuff
#include "NeugenInterface/neugen_wrapper.h"
#include "NeugenInterface/neugen_config.h"
#include "NeugenInterface/interaction.h"
#include "NeugenInterface/process.h"
#include "NeugenInterface/kinematic_variable.h"
#include "NeugenInterface/init_state.h"
#include "NeugenInterface/final_state.h"
#include "NeugenInterface/flavor.h"
#include "TLorentzVector.h"

#include "TFile.h"
#include "TROOT.h"
#include "TObject.h"
#include "TDirectory.h"
#include "TH1.h"
#include "TMath.h"
#include "TGraph.h"

#include <cstdlib>
#include <stdexcept>

CVSID("$Id: NuTransition.cxx,v 1.10 2009/10/01 19:01:24 bckhouse Exp $");

//---------------------------------------------------------------------------------
NuTransition::NuTransition(const NuXMLConfig *xmlConfig) 
{
    // Get our local xmlConfig copy, if it is safe to do so
    if (!xmlConfig) {
        // If the XML file does not exist, we can't do anything
        doNothing = true;
        
        MAXMSG("NuTransition",Msg::kInfo,5)
        << "No NuXMLConfig object supplied: "
        << "performing no systematic shifts, oscillations, or transitions."
        << endl;
        return;
    }
    
    fxmlConfig = (NuXMLConfig*)xmlConfig->Clone();
    firstRunTrans = true;

    MSG("NuTransition",Msg::kInfo) << "Applying parameters: " 
    << "dm2nu = " << xmlConfig->DM2Nu()   << ", sn2nu = " << xmlConfig->SN2Nu()
    << ", dm2bar = " << xmlConfig->DM2Bar() << ", sn2bar = " << xmlConfig->SN2Bar()
    << ", alpha = " << xmlConfig->TransitionProb() << endl;
    
    doNothing  = false;
    
    if (doNothing) return;
    
    fSyst = new NuSystematic(*fxmlConfig);
    
    // Get the binning scheme set up
    const NuUtilities cuts;
  
    //Binning scheme
    NuBinningScheme::NuBinningScheme_t binningScheme =
        static_cast<NuBinningScheme::NuBinningScheme_t>(xmlConfig->BinningScheme());
  
    std::vector<Double_t> vReco = cuts.RecoBins(binningScheme);
    std::vector<Double_t> vTrue = cuts.TrueBins(binningScheme);
    
    Int_t numRecoBins = vReco.size() - 1;
    Int_t numTrueBins = vTrue.size() - 1;
    
    Double_t *recoBinsArray = &(vReco[0]);
    Double_t *trueBinsArray = &(vTrue[0]);
    
    // Creat the reweighting histograms
    numu_flux   = new TH1D("numu_flux","numu_flux",numTrueBins,trueBinsArray);
    nubar_flux  = new TH1D("nubar_flux","nubar_flux",numRecoBins,recoBinsArray);
    flux_weight = new TH1D("flux_weight","flux_weight",numRecoBins,recoBinsArray);
    
    // Determine how the systematics are applied
    if (NuSystematic::FluxSyst(*xmlConfig)) { 
        MSG("NuTransition", Msg::kInfo) << "Not applying " << xmlConfig->Name() 
        << " to transitions b/c it's a flux systematic." << endl;
        applyTranSyst = false;
    }
    else {
        applyTranSyst = true;
    }
    
    // Get the XSec graphs from file
    TDirectory *start = gDirectory;
    
    MSG("NuTransition",Msg::kInfo) << "Get cross-section graphs from file." << endl;
    
    const char * fxsecfilename = getenv("xf");
    
    MSG("NuTransition",Msg::kInfo) << "Opening xsec file: " << fxsecfilename << endl;
    TFile *xsecfile = new TFile(fxsecfilename,"READ");
  
    // Check to see if this failed
    if (!xsecfile->IsOpen()) {
      MSG("NuTransition", Msg::kError) << "\n\n"
        << "*****************************************************************************\n"
        << "Failed to open cross section file.\n"
        << "Currently the cross section file is found through the 'xf' environment variable.\n"
        << "Perhaps you need to issue a command such as:\n   "
        << "   export xf=\"$SRT_PRIVATE_CONTEXT/NtupleUtils/data/xsec_minos_modbyrs4_v3_5_0_mk.root\"\n"
        << "*****************************************************************************\n\n"
        << endl;
      throw runtime_error("Failed to read crosssection file");
    }
    
    // NuMu
    TH1F *fhNuMuCCCrossSections = (TH1F*) xsecfile->Get("h_numu_cc_tot");
    Float_t *x = new Float_t[fhNuMuCCCrossSections->GetNbinsX()];
    Float_t *y = new Float_t[fhNuMuCCCrossSections->GetNbinsX()];
    for(int i=0;i<fhNuMuCCCrossSections->GetNbinsX();i++) {
        x[i] = fhNuMuCCCrossSections->GetBinCenter(i+1);
        y[i] = fhNuMuCCCrossSections->GetBinContent(i+1);
    }
    fNuMuXGraph = new TGraph(fhNuMuCCCrossSections->GetNbinsX(),x,y);
    if (x) {delete[] x; x = 0;}
    if (y) {delete[] y; y = 0;}
    
    // NuMuBar
    TH1F *fhNuBarCCCrossSections = (TH1F*) xsecfile->Get("h_numubar_cc_tot");
    x = new Float_t[fhNuBarCCCrossSections->GetNbinsX()];
    y = new Float_t[fhNuBarCCCrossSections->GetNbinsX()];
    for(int i=0;i<fhNuBarCCCrossSections->GetNbinsX();i++) {
        x[i] = fhNuBarCCCrossSections->GetBinCenter(i+1);
        y[i] = fhNuBarCCCrossSections->GetBinContent(i+1);
    }
    fNuBarXGraph = new TGraph(fhNuBarCCCrossSections->GetNbinsX(),x,y);
    if (x) {delete[] x; x = 0;}
    if (y) {delete[] y; y = 0;}
    
    xsecfile->Close();
    if (xsecfile){delete xsecfile; xsecfile = 0;}
    
    gDirectory = start;

}

//---------------------------------------------------------------------------------
NuTransition::~NuTransition() 
{
    if (numu_flux) delete numu_flux; numu_flux = 0;
    if (nubar_flux) delete nubar_flux; nubar_flux = 0;
    if (flux_weight) delete flux_weight; flux_weight = 0;
}

//---------------------------------------------------------------------------------
void NuTransition::Fill(NuMCEvent &mc)
{
    if (doNothing) return;
    
    // CC-only
    if (mc.iaction==1) { 
        double fw = 1;
        
        // Correct for Cross section
        fw *= this->GetXSecGraph(mc);
        
        if (mc.inu == 14) {
            numu_flux->Fill(mc.neuEnMC, mc.rw/fw);
            MAXMSG("NuTransition",Msg::kInfo,5) << "===> Fill numu, xsec=" << fw << endl;
        }
        else if (mc.inu == -14) {
            nubar_flux->Fill(mc.neuEnMC, mc.rw/fw);
            MAXMSG("NuTransition",Msg::kInfo,5) << "===> Fill nubar, xsec=" << fw << endl;
        }
        else{
            MAXMSG("NuTransition",Msg::kInfo,5) << "===> Fill nothing, ntyp=" << mc.inu << endl;
        }
    }
}

//......................................................................
void NuTransition::Reweight(NuEvent& nu)
{
    // With no xml file, do neither systematics or oscillations
    if (doNothing) return;    
 
    double wTotal = 1;
    
    // If Systematics are applied to transitioned events,
    // can simply shift the whole event
    if (applyTranSyst) {
        this->DoSystematicShifts(nu);
    }
    else {
        NuEvent * nuCopy = (NuEvent*)nu.Clone();
        this->DoSystematicShifts(*nuCopy);
        wTotal *= nuCopy->rw / nu.rw;
        delete nuCopy;
    }
    
  //bool isNuE = (nu.inu == 12 || nu.inu == -12);
    bool isNuMu = (nu.inu == 14 || nu.inu == -14);
    bool isNuTau = (nu.inu == 16 || nu.inu == -16);
  
  
    if (SimFlag::kData == nu.simFlag) {
        // don't oscillate real data
        MAXMSG("NuTransition",Msg::kInfo,1) 
        << "Not applying fake oscillations to data" << endl;
    }
    else if (Detector::kFar != nu.detector) { 
        // don't oscillate the ND data
        MAXMSG("NuTransition",Msg::kInfo,1) 
        << "Not applying fake oscillations to ND MC" << endl;
    }
    else if (1 != nu.iaction) { 
        // check for !CC (i.e. NC)
        if(isNuTau) {
            MAXMSG("NuTransition",Msg::kInfo,1) << "Removing NC events from tau file (nu.rw set to zero)" << endl;
            wTotal = 0;
        }
        MAXMSG("NuTransition",Msg::kInfo,1) << "Not oscillating NCs" << endl;
    }
    else if ( !(isNuMu || isNuTau) ) { 
        // select everything that is not NuMu(bar) or NuTau(bar) i.e. true if nue(bar)
        MAXMSG("NuTransition",Msg::kInfo,1)
        << "Not oscillating nue events (non-muon/tau (anti)neutinos)" << endl;
    }
    else if ( isNuTau && (nu.inunoosc == 12 || nu.inunoosc==-12) ) {
        // select events that started as nues and were turned in to nutaus
        MAXMSG("NuTransition",Msg::kInfo,1)
        << "Removing beam nue events from tau file (nu.rw set to zero)" << endl;
        wTotal = 0;
    }
    else if (fxmlConfig->DM2Nu() < 0.0 ||  // Do this last, so that taus and nue's are still properly
             fxmlConfig->SN2Nu() < 0.0 ||  // reweighted, even in the no-oscillation case
             fxmlConfig->DM2Bar() < 0.0 || // For example, in the no-osc case the nutaus are given a weight of 0
             fxmlConfig->SN2Bar() < 0.0){  // so have to explictly just not touch them
        MAXMSG("NuTransition",Msg::kInfo,1)
        <<"Not applying fake oscillations due to xml configuration"<<endl;
    }
    else {
        //include the oscillation weight in the general nu.rw variable
        wTotal *= this->OscillationWeight(nu.neuEnMC, nu.inu);
        MAXMSG("NuTransition",Msg::kInfo,10)
        << "Oscillating with dm2=" << fxmlConfig->DM2Nu()
        << "; weight=" << this->OscillationWeight(nu.neuEnMC, nu.inu)
        << "; energy=" << nu.neuEnMC << endl;
        
        // Only create transitioned events if there are transitions
        if (fxmlConfig->TransitionProb() > 0) {
            if (isNuTau) {
                MAXMSG("NuTransition",Msg::kInfo,1)
                << "Weighting down appeared NuTaus for transitions by " 
                << (1.0 - fxmlConfig->TransitionProb()) << endl;
                wTotal *= (1.0 - fxmlConfig->TransitionProb());
            }
            else if (nu.inu == -14) { // and this is an antineutrino
                MAXMSG("NuTransition",Msg::kInfo,10) << "Getting a transitioned event" << endl;
                
                if (firstRunTrans) {
                    flux_weight->Divide(numu_flux, nubar_flux);
                    firstRunTrans = false;
                }            
                
                // Disappearance Probability = 1 - survival probability
                Double_t transitionWeight =  1.0 - this->OscillationWeight(nu.neuEnMC,nu.inu); 
                transitionWeight *= fxmlConfig->TransitionProb();
                transitionWeight *= NuUtilities::ValueAt(flux_weight, nu.neuEnMC);
                
                // Add, not multiply, the transition weight
                wTotal += transitionWeight;
            }
        }
    }
    
    nu.rw *= wTotal;
}

/*
//---------------------------------------------------------------------------------
NuEvent* NuTransition::GetTransitionEvent(const NuEvent &nu)
{
  MAXMSG("NuTransition",Msg::kInfo,10) << "Getting a transitioned event" << endl;
  if (firstRunTrans) {
    flux_weight->Divide(numu_flux, nubar_flux);
    firstRunTrans = false;
  }
  
  
  // Does the event transition?  If not, return an empty pointer.
  if (nu.inu != -14 || nu.iaction != 1 ||
      fxmlConfig->TransitionProb() <= 0 || nu.detector != Detector::kFar){
    return 0;
  }
  
  // Make sure the flux weights have been calculated
  
  NuEvent *nuTrans = (NuEvent*)nu.Clone();
  
  if (applyTranSyst) {
    this->DoSystematicShifts(*nuTrans);
  }
  
  // Disappearance Probability = 1 - survival probability
  Double_t transitionWeight =  1.0 - this->OscillationWeight(nuTrans->neuEnMC,nuTrans->inu,fxmlConfig);
  transitionWeight *= fxmlConfig->TransitionProb();
  transitionWeight *= NuUtilities::ValueAt(flux_weight, nuTrans->neuEnMC);
  
  nuTrans->rw *= (transitionWeight);
  
  if (nuTrans->rw == 0) MSG("NuTransition",Msg::kWarning) << "Warning: zero weight" << endl;
  
  return nuTrans;
}           
*/

//---------------------------------------------------------------------------------
void NuTransition::DoSystematicShifts(NuEvent& event)
{
    MAXMSG("NuTransition",Msg::kInfo,5)
    << "Performing systematic shift "
    << fxmlConfig->FullTitle()
    << endl;
    fSyst->Shift(event);
}

//......................................................................

const Float_t NuTransition::OscillationWeight
(const Float_t energy, const Int_t inu, const NuXMLConfig* lxmlConfig) const
{
    const NuXMLConfig *xmlConfig;
    if (lxmlConfig) xmlConfig = lxmlConfig;
    else            xmlConfig = fxmlConfig;
        
    if (14==inu){
      return NuUtilities::OscillationWeight(energy,
                                            xmlConfig->DM2Nu(),
                                            xmlConfig->SN2Nu());
    }
    else if (-14==inu){
      return NuUtilities::OscillationWeight(energy,
                                            xmlConfig->DM2Bar(),
                                            xmlConfig->SN2Bar());
    }
    if (16==inu){
      return 1-NuUtilities::OscillationWeight(energy,
                                              xmlConfig->DM2Nu(),
                                              xmlConfig->SN2Nu());
    }
    else if (-16==inu){
      return 1-NuUtilities::OscillationWeight(energy,
                                              xmlConfig->DM2Bar(),
                                              xmlConfig->SN2Bar());
    }
    else return 1; // Doesn't oscillate
}

//......................................................................

Double_t NuTransition::GetXSecGraph(NuEvent& nu, bool opp_charge) const
{
    Double_t nuBarXSec = fNuBarXGraph->Eval(nu.neuEnMC,0,"");
    Double_t nuMuXSec =  fNuMuXGraph->Eval( nu.neuEnMC,0,"");
    
    bool numu = (nu.inu > 0);
    if (opp_charge) numu = !numu;
    
    if (numu) return nuMuXSec;
    else      return nuBarXSec;
}
//......................................................................

Double_t NuTransition::GetXSecGraph(NuMCEvent& nu, bool opp_charge) const
{
    Double_t nuBarXSec = fNuBarXGraph->Eval(nu.neuEnMC,0,"");
    Double_t nuMuXSec =  fNuMuXGraph->Eval( nu.neuEnMC,0,"");
    
    bool numu = (nu.inu > 0);
    if (opp_charge) numu = !numu;
       
    if (numu) return nuMuXSec;
    else      return nuBarXSec;
}

//......................................................................

Double_t NuTransition::GetXSecNeugen(NuEvent& nu, bool opp_charge) const
{
    // Get Cross Sections //
    
    static neugen_wrapper *fNuWrap = 0;
    //static Registry std_registry;   
    static neugen_config *std_config = 0;
    
    if (!fNuWrap) {
        //std_registry.Set("neugen:config_name","MODBYRS");
        //std_registry.Set("neugen:config_no",4);
        
        //NuSystematic::SetNeugenDefaults(std_registry);
        
        std_config = new neugen_config("stdNeuConfig");
        std_config->set_best_parameters();
        //if(config) SetStandardConfig(config);
        
        fNuWrap = new neugen_wrapper(std_config);
    }
    
    //event kinematic variables:
    double kval[5] = {0,0,0,0,0};
    int kid1 = -1;
    int kid2 = -1;
    kval[1] = nu.q2MC;
    kval[2] = nu.w2MC; 
    kval[3] = nu.xMC;
    kval[4] = nu.yMC;
    
    //translate from minossoft-speak to neugen-speak:
    int flavor = 4;
    if(abs(nu.inu)==12) flavor = 1;
    else if(abs(nu.inu)==14) flavor = 2; 
    else if(abs(nu.inu)==16) flavor = 3;
    else if(abs(nu.inu)==1 || abs(nu.inu)==2 || abs(nu.inu)==3) flavor = abs(nu.inu);
    else {
        MSG("NuTransition",Msg::kError) << "Invalid event:inu value " << nu.inu << endl;
        return 0;
    }
    
    int ccnc = nu.iaction;
    if(ccnc==0) ccnc = 2;
    if(!(ccnc==1||ccnc==2)) {
        MSG("NuTransition",Msg::kError) << "Invalid event:iaction value " 
        << nu.iaction << endl;
        return 0;
    }
    
    int process = 0;
    if(nu.iresonance>=1001&&nu.iresonance<=1004) process = nu.iresonance - 1000;
    else if(nu.iresonance>=1&&nu.iresonance<=4) process = nu.iresonance;
    else {
        if(nu.iresonance==1005) MSG("NuTransition",Msg::kInfo) 
            << "Event:iresonance value is 1005 "
            << "- no reweighting performed" << endl;
        else MSG("NeugenWC",Msg::kError) 
            << "Invalid event:iresonance (or event:process) value "
            << nu.iresonance << endl;
        return 0;
    }
    
    //extract necessary kinematic variables:
    if(process==1) {   //qel
        kid1=1;
        kid2=0;
        if(kval[1]==999999) {
            MSG("NuTransition",Msg::kWarning) 
            << "Current event is QEL: q^2 kinematic variable not set!"
            << " Returning weight of 1." << endl;
            return 0;
        }    
    }
    else if(process==2) { //res
        kid1=1;
        kid2=2;
        if(kval[1]==999999) {
            MSG("NuTransition",Msg::kWarning) 
            << "Current event is RES: q^2 kinematic variable not set!"
            << " Returning weight of 1." << endl;
            return 0;
        }
        if(kval[2]==999999) {
            MSG("NuTransition",Msg::kWarning) 
            << "Current event is RES: w^2 kinematic variable not set!"
            << " Returning weight of 1." << endl;
            return 0;
        }
        //need to pass w rather than w^2 to neugen
        if(kval[2]>0) kval[2]=sqrt(kval[2]); 
        else {
            MSG("NuTransition",Msg::kWarning) 
            << "w^2 kinematic variable is negative!"
            << " Returning weight of 1." << endl;
            return 0;
        }
    }
    else if(process==3) { //dis
        kid1=3;
        kid2=4;
        if(kval[3]==999999) {
            MSG("NuTransition",Msg::kWarning) 
            << "Current event is DIS: x kinematic variable not set!"
            << " Returning weight of 1." << endl;
            return 0;
        }
        if(kval[4]==999999) {
            MSG("NuTransition",Msg::kWarning) 
            << "Current event is DIS: y kinematic variable not set!"
            << " Returning weight of 1." << endl;
            return 0;
        }
    }
    else if(process==4) { //coh
        //reweighting not possible here yet
        return 0;
    }
    
    //start making neugen style objects:
    TLorentzVector *nu_lv = new TLorentzVector(nu.neuPxMC,nu.neuPyMC,nu.neuPzMC,nu.neuEnMC);
    TLorentzVector *tar_lv = new TLorentzVector(nu.tgtPxMC,nu.tgtPyMC,nu.tgtPzMC,nu.tgtEnMC);
    
    init_state_t init = init_state_enum(nu.initialStateMC);
    
    if (opp_charge) {
        using namespace init_state;
        switch(init) {
            case e_vN:   init = e_vbN;  break;
            case e_vbN:  init = e_vN;   break;
                
            case e_vA:   init = e_vbA;  break;
            case e_vbA:  init = e_vA;   break;
                
            case e_vp:   init = e_vbp;  break;
            case e_vbp:  init = e_vp;   break;
                
            case e_vn:   init = e_vbn;  break;
            case e_vbn:  init = e_vn;   break;
                
            case e_undefined_init_state:
            default:            
                MSG("NuTransition",Msg::kWarning) << "Unknown initial state" << endl;
                return 0;
        }
        
        if (process == 1) {
            if (init == e_vbn) init = e_vbp;
            if (init == e_vp)  init = e_vn;
        }  
    }
    
    interaction *iact = new interaction(flavor_enum(flavor),
                                        nucleus_enum(nu.nucleusMC),
                                        ccnc_enum(ccnc),
                                        init); 
    
    final_state *fs = new final_state(nu.hadronicFinalStateMC);
    
    
    //calculate cross section
    Double_t sigma = fNuWrap->offshell_diff_xsec(nu_lv,tar_lv,
                                                 iact,process_enum(process),fs,
                                                 kinematic_variable_enum(kid1),float(kval[kid1]),
                                                 kinematic_variable_enum(kid2),float(kval[kid2]));
    
    
    delete nu_lv;
    delete tar_lv;
    delete iact;
    delete fs;
    
    
    return sigma;
}

---