This notebook calculates the chi2 for the MINERvA QElike data release directly from the website.
Models for qelike are:

Models={"MnvGENIE","MnvGENIE_RPA_2p2h_notune","MnvGENIE_RPA_no2p2h","MnvGENIE_noRPA_2p2h","MnvGENIE_noRPA_no2p2h","Untuned_GENIE","gfg_norpa_Nieves_ma099","gfg_norpa_nomec_ma099","gfg_norpa_tem_ma099","gfg_rpa_Nieves_ma099","gfg_rpa_tem_ma099","lfg_norpa_tem_ma099","mc_sf_norpa_nomec_ma099"};

model = Models[[2]]

Units for cross sections are 10^-41 cm^2

SETUP

Get the data and MC - change directory to location of csv files

signal = "qelike";

source = StringJoin["http://physics.oregonstate.edu/~schellmh/",signal,"/"];

header =StringJoin[ source,"cross_sections_muonpz_muonpt_lowangleqelike_minerva"];
mcheader = StringJoin[source,"MnvQElike-pzpt-Model-"];
tail := "_2d.csv";
If[signal == "ccqe",{header= StringJoin[ source,"cross_sections_muonpz_muonpt_lowangleccqe2p2h_minerva"],
mcheader = StringJoin[source,"MnvCCQE-pzpt-Model-"]}
];

MyPlot[array_,range_]:=ListContourPlot[array,PlotLegends->Automatic,InterpolationOrder->0,ImageSize->300,PlotRange->range]

Get the files

centralname = StringJoin[header,"",tail]

central = Import[centralname]

mcname = StringJoin[mcheader,model,tail]

mc = Import[mcname]

statname = StringJoin[header,"_staterrors",tail]

staterrs = Import[statname]

sysname = StringJoin[header,"_syserrors",tail]

syserrs = Import[sysname]

covname = StringJoin[header,"_covariance.csv"]

thecovariance =  Import[covname];

Need to flatten the 2 d info to use the covariance matrix

centraldata =  Flatten[central]

and eliminate the bins with zero entries

goodelements = Flatten[SparseArray[centraldata]["NonzeroPositions"]]

len = Length[goodelements]

central1d = Table[centraldata[[goodelements[[i]]]],{i,1,len}]

mc1d = Table[Flatten[mc][[goodelements[[i]]]],{i,1,len}]

residual = central1d - mc1d

staterrs1d = Table[Flatten[staterrs][[goodelements[[i]]]],{i,1,len}];

syserrs1d = Table[Flatten[syserrs][[goodelements[[i]]]],{i,1,len}];

toterr = Sqrt[staterrs1d^2 + syserrs1d^2];

covariance= Table[thecovariance[[goodelements[[i]],goodelements[[j]]]],{i,1,len},{j,1,len}] ;

Do some checks for roundoff error

diag =Sqrt[Diagonal[covariance]]

toterr/diag

correlation = Table[covariance[[i,j ]]/diag[[i]]/diag[[j]],{i,1,len},{j,1,len}];

MyPlot[correlation,{-1,1}]

Linear Chi2

InvCovariance = Inverse[covariance];

MyPlot[InvCovariance,Full]

chi2cont = Table[InvCovariance[[i]][[j]]*residual[[i]]*residual[[j]],{i,1,Length[residual]},{j,1,Length[residual]}];

MyPlot[chi2cont,Full]

chi2 = Sum[InvCovariance[[i]][[j]]*residual[[i]]*residual[[j]],{i,1,Length[residual]},{j,1,Length[residual]}]

Log Version

Make the logs

LogCentral = Log[central1d]

LogMC = Log[mc1d]

LogResidual =LogCentral - LogMC

LogCovariance = Table[covariance[[i,j]]/central1d[[i]]/central1d[[j]],{i,1,len},{j,1,len}];

LogInvCovariance = Inverse[LogCovariance];

LogChi2cont = Table[LogInvCovariance[[i]][[j]]*LogResidual[[i]]*LogResidual[[j]],{i,1,len},{j,1,len}];

MyPlot[LogChi2cont,Full]

LogChi2 = Sum[LogInvCovariance[[i]][[j]]*LogResidual[[i]]*LogResidual[[j]],{i,1,len},{j,1,len}]

Result - approx VI

{model,chi2,LogChi2}