diff --git a/data/flux/ANL_1977_2horn_rescan.root b/data/flux/ANL_1977_2horn_rescan.root
new file mode 100644
index 0000000..3524d3a
Binary files /dev/null and b/data/flux/ANL_1977_2horn_rescan.root differ
diff --git a/data/flux/BNL_NuInt02_rescan.root b/data/flux/BNL_NuInt02_rescan.root
new file mode 100644
index 0000000..ed4604e
Binary files /dev/null and b/data/flux/BNL_NuInt02_rescan.root differ
diff --git a/parameters/config.xml b/parameters/config.xml
index 834932e..6507c32 100644
--- a/parameters/config.xml
+++ b/parameters/config.xml
@@ -1,187 +1,192 @@
<nuisance>
<!-- # ###################################################### -->
<!-- # NUISANCE CONFIGURATION OPTIONS -->
<!-- # This file is read in by default at runtime -->
<!-- # If you want to override on a case by case bases use -q at runtime -->
<!-- # ###################################################### -->
<!-- # MAIN Configs -->
<!-- # ###################################################### -->
<!-- # Logger goes from -->
<!-- # 1 Quiet -->
<!-- # 2 Fitter -->
<!-- # 3 Samples -->
<!-- # 4 Reconfigure Loops -->
<!-- # 5 Every Event print out (SHOUT) -->
<!-- # -1 DEBUGGING -->
<config verbosity='5'/>
<config VERBOSITY='5'/>
<!-- # ERROR goes from -->
<!-- # 0 NONE -->
<!-- # 1 FATAL -->
<!-- # 2 WARN -->
<config ERROR='2'/>
<config TRACE='0'/>
<config cores='2' />
<config spline_test_throws='50' />
<config spline_cores='6' />
<config spline_chunks='10' />
<config spline_procchunk='-1' />
<config Electron_NThetaBins='4' />
<config Electron_NEnergyBins='4' />
<config Electron_ThetaWidth='1.0' />
<config Electron_EnergyWidth='0.10' />
<config RemoveFSIParticles='0' />
<config RemoveUndefParticles='0' />
<config RemoveNuclearParticles='0'/>
<config logging.JointFCN.cxx='4'/>
<!-- # Input Configs -->
<!-- # ###################################################### -->
<!-- # Default Requirements file for the externalDataFitter Package -->
<!-- # MAX Events : -1 is no cut. Events will be scaled automatically to give good xsec predictions. -->
<config input.maxevents='-1'/>
<config MAXEVENTS='-1'/>
<config input.MAXEVENTS='-1'/>
<config includeemptystackhists='0'/>
<!-- # Turn on/off event manager -->
<!-- # EventManager enables us to only loop number of events once for multiple projections of the same measurements -->
<!-- # e.g. MiniBooNE CC1pi+ Q2 and MiniBooNE CC1pi+ Tmu would ordinarily require 2 reconfigures, but with this enabled it requires only one -->
<config input.eventmanager='1'/>
<config EventManager='1'/>
<!-- # Event Directories -->
<!-- # Can setup default directories and use @EVENT_DIR/path to link to it -->
<config EVENT_DIR='/data2/stowell/NIWG/'/>
<config NEUT_EVENT_DIR='/data2/stowell/NIWG/neut/fit_samples_neut5.3.3/'/>
<config GENIE_EVENT_DIR='/data2/stowell/NIWG/genie/fit_samples_R.2.10.0/'/>
<config NUWRO_EVENT_DIR='/data2/stowell/NIWG/nuwro/fit_samples/'/>
<config GIBUU_EVENT_DIR='/data/GIBUU/DIR/'/>
<config SaveNuWroExtra='0' />
<!-- # In PrepareGENIE the reconstructed splines can be saved into the file -->
<config save_genie_splines='1'/>
<!-- # In InputHandler the option to regenerate NuWro flux/xsec plots is available -->
<!-- # Going to move this to its own app soon -->
<config input.regen_nuwro_plots='0'/>
<!-- # DEVEL CONFIG OPTION, don't touch! -->
<config CacheSize='0'/>
<!-- # ReWeighting Configuration Options -->
<!-- # ###################################################### -->
<!-- # Set absolute twkdial for parameters -->
<config params.setabstwk='0'/>
<!-- # Convert Dials in output statements using dial conversion card -->
<config convert_dials='0'/>
<!-- # Make RW Calculations be quiet -->
<config params.silentweighting='0'/>
<!-- # Vetos can be used to specify RW dials NOT to be loaded in -->
<!-- # Useful if one specific one has an issue -->
<config FitWeight.fNIWGRW_veto=''/>
<config FitWeight.fNuwroRW_veto=''/>
<config FitWeight.fNeutRW_veto=''/>
<config FitWeight.fGenieRW_veto=''/>
<!-- # Output Options -->
<!-- # ###################################################### -->
<!-- # Save Nominal prediction with all rw engines at default -->
<config savenominal='0'/>
<!-- # Save prefit with values at starting values -->
<config saveprefit='0'/>
<!-- # Here's the full list of drawing options -->
<!-- # See src/FitBase/Measurement1D::Write for more info -->
<!-- #config drawopts DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/RESIDUAL/MATRIX/FLUX/MASK/MAP -->
<!-- #config drawopts DATA/MC -->
<config drawopts='DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/FLUX/XSEC/MASK/COV/INCOV/DECOMP/CANVPDG/CANVMC/SETTINGS'/>
<config InterpolateSigmaQ0Histogram='1' />
<config InterpolateSigmaQ0HistogramRes='100' />
<config InterpolateSigmaQ0HistogramThrow='1' />
<config InterpolateSigmaQ0HistogramNTHROWS='100000' />
<!-- # Save the shape scaling applied with option SHAPE into the main MC hist -->
<config saveshapescaling='0'/>
<config CorrectGENIEMECNorm='1'/>
<!-- # Set style of 1D output histograms -->
<config linecolour='1'/>
<config linestyle='1'/>
<config linewidth='1'/>
<!-- # For GenericFlux -->
<config isLiteMode='0'/>
<!-- # Statistical Options -->
<!-- # ###################################################### -->
<!-- # Add MC Statistical error to likelihoods -->
<config statutils.addmcerror='0'/>
<!-- # NUISMIN Configurations -->
<!-- # ###################################################### -->
<config minimizer.maxcalls='1000000'/>
<config minimizer.maxiterations='1000000'/>
<config minimizer.tolerance='0.001'/>
<!-- # Number of events required in low stats routines -->
<config minimizer.lowstatevents='25000'/>
<!-- # Error band generator configs -->
<!-- # ###################################################### -->
<!-- # For -f ErrorBands creates error bands for given measurements -->
<!-- # How many throws do we want (The higher the better precision) -->
<config error_throws='250'/>
<!-- # Are we throwing uniform or according to Gaussian? -->
<!-- # Only use uniform if wanting to study the limits of a dial. -->
<config error_uniform='0'/>
<config WriteSeperateStacks='1'/>
<!-- # Other Individual Case Configs -->
<!-- # ###################################################### -->
<!-- # Covariance throw options for fake data studies with MiniBooNE data. -->
<config thrown_covariance='FULL'/>
<config throw_mc_stat='0.0'/>
<config throw_diag_syst='0'/>
<config throw_corr_syst='0'/>
<config throw_mc_stat='0'/>
<!-- # Apply a shift to the muon momentum before calculation of Q2 -->
<config muon_momentum_shift='0.0'/>
<config muon_momentum_throw='0'/>
<!-- # MINERvA Specific Configs -->
<config MINERvA_XSec_CCinc_2DEavq3_nu.hadron_cut='0'/>
<config MINERvA_CCinc_XSec_2DEavq3_nu.useq3true='0'/>
<config Modes.split_PN_NN='0'/>
<config SignalReconfigures='0'/>
<!-- # SciBooNE specific -->
<config SciBarDensity='1.04'/>
<config SciBarRecoDist='10.0'/>
<config PenetratingMuonEnergy='1.4'/>
<config NumRangeSteps='50'/>
+
+<config GENIEWeightEngine_CCRESMode="kModeMaMv"/>
+<config GENIEInputHandler.SavePrimary="0" />
+
</nuisance>
+
diff --git a/scripts/plotlimitdist.py b/scripts/plotlimitdist.py
new file mode 100644
index 0000000..cbb6134
--- /dev/null
+++ b/scripts/plotlimitdist.py
@@ -0,0 +1,178 @@
+from ROOT import *
+from array import array
+import sys
+
+def GetUID(name):
+ return name.replace("_" + name.split("_")[-1],"")
+def GetValue(name):
+ return float(name.split("_")[-1])
+
+
+colorlist = [kGreen,kRed, kBlue, kOrange]
+
+infile = TFile(sys.argv[1], "READ")
+alllimitdir = infile.Get("Limits")
+
+# Make starting legend
+leg = TLegend(0.1,0.1,0.9,0.9)
+alluidstyles = {}
+alllines = []
+
+c1 = TCanvas("c1","c1",1000,600)
+
+for dirkey in alllimitdir.GetListOfKeys():
+ if dirkey.GetName() == "nominal":
+ linedata = TLine(0.1,0.1,0.9,0.9)
+ linedata.SetLineColor(kBlack)
+ leg.AddEntry(linedata, "DATA", "l")
+ alllines.append(linedata)
+
+ linemc = TLine(0.1,0.1,0.9,0.9)
+ linemc.SetLineColor(kGreen)
+ leg.AddEntry(linemc, "Nominal MC", "l")
+ alllines.append(linemc)
+ alluidstyles["nominal"] = kGreen
+
+ else:
+ uid = GetUID(dirkey.GetName())
+ if uid not in alluidstyles:
+ alluidstyles[uid] = colorlist[len(alluidstyles)]
+ line = TLine(0.1,0.1,0.9,0.9)
+ line.SetLineColor( alluidstyles[uid] )
+ leg.AddEntry(line, uid, "l")
+ alllines.append(line)
+
+leg.Draw()
+gPad.Update()
+c1.SaveAs("limitplots.pdf(")
+gStyle.SetOptTitle(1)
+c1.Clear()
+c1.Divide(3,1)
+
+# Get Nominal
+nomdir = infile.Get("Limits/nominal")
+for plotkey in nomdir.GetListOfKeys():
+
+ # Only get data
+ if (not plotkey.GetName().endswith("_data")): continue
+
+ # Get Plots
+ data = nomdir.Get(plotkey.GetName())
+ datashape = nomdir.Get(plotkey.GetName()).Clone()
+ mc = nomdir.Get(plotkey.GetName().replace("data","MC"))
+ mcshape = nomdir.Get(plotkey.GetName().replace("data","MC_SHAPE"))
+
+ if not data or not mc: continue
+
+ # Draw nominal and data
+ c1.cd(1)
+ data.GetYaxis().SetRangeUser(0.0, data.GetMaximum()*1.5)
+ data.SetLineColor(kBlack)
+ data.SetLineWidth(3)
+ data.Draw("E1")
+ mc.SetLineColor(kGreen)
+ mc.SetLineWidth(3)
+ mc.Draw("SAME HIST")
+ alllimitmc = []
+
+ c1.cd(2)
+ datashape.GetYaxis().SetRangeUser(0.0, data.GetMaximum()*1.5)
+ datashape.SetLineColor(kBlack)
+ datashape.SetLineWidth(3)
+ datashape.SetTitle(datashape.GetName() + "_SHAPE")
+ datashape.Draw("E1")
+ mcshape.SetLineWidth(3)
+ mcshape.SetLineColor(kGreen)
+ mcshape.Draw("SAME HIST")
+
+ likelihoodvals = {}
+ maxlike = 0.0
+ minlike = 1.E9
+
+
+
+ # Loop through and get limit plots
+ for dirkey in alllimitdir.GetListOfKeys():
+ if dirkey.GetName() == "nominal": continue
+ limitdir = alllimitdir.Get(dirkey.GetName())
+ limitmc = limitdir.Get(plotkey.GetName().replace("_data","_MC")).Clone()
+ limitmcshape = limitdir.Get(plotkey.GetName().replace("_data","_MC_SHAPE")).Clone()
+
+ print "Getting from key:", limitdir, limitmc
+ uid = GetUID(dirkey.GetName())
+ if uid not in likelihoodvals:
+ likelihoodvals[uid] = []
+
+ limitmc.SetLineColor( alluidstyles[uid] )
+ limitmc.SetLineWidth(1)
+
+ limitmcshape.SetLineColor( alluidstyles[uid] )
+ limitmcshape.SetLineWidth(1)
+
+ c1.cd(1)
+ limitmc.Draw("SAME HIST C")
+ alllimitmc.append(limitmc)
+
+ c1.cd(2)
+ limitmcshape.Draw("SAME HIST C")
+ alllimitmc.append(limitmcshape)
+
+ xval = float(GetValue(dirkey.GetName()))
+ yval = float(limitmc.GetTitle())
+ likelihoodvals[uid].append([xval,yval])
+
+ maxlike = max([maxlike, (float(limitmc.GetTitle()))])
+ minlike = min([minlike, (float(limitmc.GetTitle()))])
+
+
+
+ #likelihoodvals["nominal"] = [[-1000, minlike],
+# [1000, minlike]]
+# alluidstyles["nominal"] = kGreen
+# likelihoodvals["nominal2"] = [[1000, minlike+1.00],
+# [-1000, minlike+1.00]]
+# alluidstyles["nominal2"] = kGreen
+# likelihoodvals["nominal3"] = [[-1000, minlike+7.82],
+# [1000, minlike+7.82]]
+# alluidstyles["nominal3"] = kGreen
+
+ # Make like hist
+ c1.cd(3)
+ dif = (maxlike - minlike)/10.0
+ allgr = []
+ for i, uid in enumerate(sorted(likelihoodvals)):
+ print i, uid
+ xvals = []
+ yvals = []
+ for valset in sorted(likelihoodvals[uid]):
+ xvals.append(valset[0])
+ yvals.append(valset[1])
+
+
+ gr = TGraph(len(xvals), array('f',xvals), array('f',yvals))
+ gr.SetTitle("Likelihood Scan about Nominal")
+ gr.SetLineColor( alluidstyles[uid] )
+
+ if "nominal" in uid:
+ gr.SetMarkerStyle(20)
+ gr.SetMarkerColor(kGreen)
+ if i == 0:
+ gr.Draw("APL")
+ else: gr.Draw("SAME PL")
+
+ gr.GetYaxis().SetRangeUser(minlike - dif, maxlike + dif)
+ allgr.append(gr)
+
+ c1.Update()
+ c1.SaveAs("limitplots.pdf")
+
+c1.SaveAs("limitplots.pdf)")
+
+
+
+
+
+
+
+
+
diff --git a/scripts/plotnuiscomp.py b/scripts/plotnuiscomp.py
new file mode 100644
index 0000000..13369a1
--- /dev/null
+++ b/scripts/plotnuiscomp.py
@@ -0,0 +1,175 @@
+from ROOT import *
+import os
+import sys
+
+
+gColorList = [kRed, kGreen, kBlue, kYellow, kOrange, kBlack]
+
+def DrawDataMC(keyname, filelist):
+
+ # Extract Data
+ data = None
+ for readfile in filelist:
+ print keyname
+ data = readfile[0].Get(keyname)
+ if not data: continue
+ break
+
+ if not data:
+ print "Data not found for : ", keyname
+ sys.exit(-1)
+
+ # Main Data Formatting
+ data.SetTitle(keyname)
+ data.SetLineColor(kBlack)
+ data.SetLineWidth(2)
+
+ # Extract MC
+ singlemclist = []
+ singledatalist = []
+ for i, mcfile in enumerate(allfiles):
+
+ print mcfile[0]
+ # Extract individual MC
+ mckey = keyname.replace("_data","_MC")
+ singlemc = mcfile[0].Get(mckey)
+ if singlemc:
+
+ singlemc = singlemc.Clone(mcfile[1]+"_MC")
+ singlemc.SetLineColor( gColorList[i] )
+ singlemc.SetLineWidth(2)
+ singlemc.SetTitle( mcfile[1] + " (" + str(singlemc.GetTitle().strip()) + ") " )
+
+ singlemclist.append(singlemc.Clone())
+ del singlemc
+
+ # Extra individual data (optional)
+ singledata = mcfile[0].Get(keyname)
+ if singledata:
+
+ singledata = singledata.Clone(mcfile[1] + "_DATA")
+ singledata.SetLineColor( kBlack )
+ singledata.SetLineWidth(2)
+ singledata.SetTitle( "^-> Saved Data" )
+
+ singledatalist.append(singledata.Clone())
+ del singledata
+
+
+ # Assign Ranges
+ miny = 99999.9
+ maxy = 0.0
+ for i in range(data.GetNbinsX()):
+ miny = min([data.GetBinContent(i+1) - data.GetBinError(i+1),miny])
+ maxy = max([data.GetBinContent(i+1) + data.GetBinError(i+1),maxy])
+ for singlemc in singlemclist:
+ miny = min([singlemc.GetMinimum(),miny])
+ maxy = max([singlemc.GetMaximum(),maxy])
+ for singledata in singledatalist:
+ miny = min([singledata.GetMinimum(),miny])
+ maxy = max([singledata.GetMaximum(),maxy])
+ widthy = maxy - miny
+
+ # Assign Ranges to data
+ if "1D" in keyname: data.GetYaxis().SetRangeUser( miny - 0.1*widthy, maxy + 0.3*widthy)
+ elif "2D" in keyname: data.GetZaxis().SetRangeUser( miny - 0.1*widthy, maxy + 0.3*widthy)
+
+ # Draw Plots 1D
+ if "1D" in keyname:
+ data.Draw("E1")
+ for mc in singlemclist:
+ mc.Draw("SAME HIST")
+
+ # Draw Plots 2D
+ elif "2D" in keyname:
+ data.Draw("E1")
+ for mc in singlemclist:
+ mc.Draw("SAME LEGO")
+
+ # Build Legend
+ leg = gPad.BuildLegend(0.45,0.65,0.8,0.85)
+ leg.SetFillStyle(0)
+ leg.SetFillColorAlpha(0,0.0)
+ leg.SetBorderSize(0)
+
+ gStyle.SetOptTitle(1)
+ gPad.SetGridx(0)
+ gPad.SetGridy(0)
+ gPad.Update()
+
+ singlemclist.append(data)
+ return singlemclist
+
+
+if __name__=="__main__":
+ c1 = TCanvas("c1","c1",210*4,297*4)
+ c1.cd()
+
+
+ # Make filelist
+ allfiles = []
+
+ for i in xrange(2, len(sys.argv)):
+ print "Reading ", i, sys.argv[i]
+
+ # split by comma
+ splitname = sys.argv[i].split(",")
+
+ # Get First
+ if (os.path.isfile(splitname[0])):
+
+ # Get File
+ newfile = (TFile(splitname[0],"READ"))
+ if not newfile:
+ print "File is not a ROOT file : ", splitname[0]
+ sys.exit()
+
+ # Set Name
+ name = splitname[0].replace(".root","")
+ if len(splitname) > 1:
+ name = splitname[1]
+
+ allfiles.append([newfile, name])
+
+
+ print allfiles
+
+ # Parse Unique Keys
+ uniquekeys = []
+ for readfile in allfiles:
+ for readkey in readfile[0].GetListOfKeys():
+ if not (readkey.GetName().endswith("_data")): continue
+ if readkey.GetName() not in uniquekeys: uniquekeys.append(readkey.GetName())
+
+ print uniquekeys
+
+ # Setup First Page
+ leg = TLegend(0.1,0.1,0.9,0.9)
+ for i, readfile in enumerate(allfiles):
+ hist = TH1D(readfile[1],readfile[1],1,0,1)
+ hist.SetLineColor(gColorList[i % len(gColorList)])
+ hist.SetLineWidth(2)
+
+ leg.AddEntry(hist, readfile[1], "l")
+
+ leg.Draw()
+ gPad.Update()
+
+
+ outputfile = sys.argv[1]
+ c1.Print(outputfile + "(")
+
+
+ # Loop through unique keys
+ for readkey in uniquekeys:
+
+ # Draw
+ datamclist = DrawDataMC(readkey, allfiles)
+
+ # Save
+ c1.Print(outputfile)
+
+ # Now save the legend again to close...
+ leg.Draw()
+ gPad.Update()
+ gPad.Print(outputfile + ")")
diff --git a/scripts/plotnuismin.py b/scripts/plotnuismin.py
new file mode 100644
index 0000000..1ef01c7
--- /dev/null
+++ b/scripts/plotnuismin.py
@@ -0,0 +1,273 @@
+from ROOT import *
+import os
+import sys
+
+
+gColorList = [kRed, kGreen, kBlue, kYellow, kOrange, kBlack]
+
+def DrawDataMC(keyname, filelist):
+
+ # Extract Data
+ data = None
+ for readfile in filelist:
+ print keyname
+ data = readfile[0].Get(keyname)
+ if not data: continue
+ break
+
+ if not data:
+ print "Data not found for : ", keyname
+ sys.exit(-1)
+
+ # Main Data Formatting
+ data.SetTitle(keyname)
+ data.SetLineColor(kBlack)
+ data.SetLineWidth(2)
+
+ # Extract MC
+ singlemclist = []
+ singledatalist = []
+ for i, mcfile in enumerate(allfiles):
+
+ print mcfile[0]
+ # Extract individual MC
+ mckey = keyname.replace("_data","_MC")
+ singlemc = mcfile[0].Get(mckey)
+ if singlemc:
+
+ singlemc = singlemc.Clone(mcfile[1]+"_MC")
+ singlemc.SetLineColor( gColorList[i] )
+ singlemc.SetLineWidth(2)
+ singlemc.SetTitle( mcfile[1] + " (" + str(singlemc.GetTitle().strip()) + ") " )
+
+ singlemclist.append(singlemc.Clone())
+ del singlemc
+
+ # Extra individual data (optional)
+ singledata = mcfile[0].Get(keyname)
+ if singledata:
+
+ singledata = singledata.Clone(mcfile[1] + "_DATA")
+ singledata.SetLineColor( kBlack )
+ singledata.SetLineWidth(2)
+ singledata.SetTitle( "^-> Saved Data" )
+
+ singledatalist.append(singledata.Clone())
+ del singledata
+
+
+ # Assign Ranges
+ miny = 99999.9
+ maxy = 0.0
+ for i in range(data.GetNbinsX()):
+ miny = min([data.GetBinContent(i+1) - data.GetBinError(i+1),miny])
+ maxy = max([data.GetBinContent(i+1) + data.GetBinError(i+1),maxy])
+ for singlemc in singlemclist:
+ miny = min([singlemc.GetMinimum(),miny])
+ maxy = max([singlemc.GetMaximum(),maxy])
+ for singledata in singledatalist:
+ miny = min([singledata.GetMinimum(),miny])
+ maxy = max([singledata.GetMaximum(),maxy])
+ widthy = maxy - miny
+
+ # Assign Ranges to data
+ if "1D" in keyname: data.GetYaxis().SetRangeUser( miny - 0.1*widthy, maxy + 0.3*widthy)
+ elif "2D" in keyname: data.GetZaxis().SetRangeUser( miny - 0.1*widthy, maxy + 0.3*widthy)
+
+ # Draw Plots 1D
+ if "1D" in keyname:
+ data.Draw("E1")
+ for mc in singlemclist:
+ mc.Draw("SAME HIST")
+
+ # Draw Plots 2D
+ elif "2D" in keyname:
+ data.Draw("E1")
+ for mc in singlemclist:
+ mc.Draw("SAME LEGO")
+
+ # Build Legend
+ leg = gPad.BuildLegend(0.45,0.65,0.8,0.85)
+ leg.SetFillStyle(0)
+ leg.SetFillColorAlpha(0,0.0)
+ leg.SetBorderSize(0)
+
+ gStyle.SetOptTitle(1)
+ gPad.SetGridx(0)
+ gPad.SetGridy(0)
+ gPad.Update()
+
+ singlemclist.append(data)
+ return singlemclist
+
+def DrawFitDialsPlot(allfiles):
+
+ singlemclist = []
+ singlelimitlist = []
+ for i, mcfile in enumerate(allfiles):
+
+ singlemc = mcfile[0].Get("fit_dials")
+ if not singlemc: continue
+
+ # Setup fit result
+ singlemc = singlemc.Clone(mcfile[1]+"_FIT")
+ singlemc.SetLineColor( gColorList[i] )
+ singlemc.SetFillColorAlpha( gColorList[i], 0.3 )
+ singlemc.SetLineWidth(2)
+ singlemc.SetTitle( mcfile[1] )
+
+ singlemclist.append(singlemc.Clone())
+ del singlemc
+
+ # Setup Limits
+ singlestart = mcfile[0].Get("start_dials")
+ singlemin = mcfile[0].Get("min_dials")
+ singlemax = mcfile[0].Get("max_dials")
+
+ print singlestart, singlemin, singlemax
+
+ singlestart.SetLineColor(gColorList[i])
+ singlestart.SetLineWidth(1)
+ singlestart.SetLineStyle(7)
+ singlelimitlist.append(singlestart.Clone())
+
+ singlemin.SetLineColor(gColorList[i])
+ singlemin.SetLineWidth(2)
+ singlemin.SetLineStyle(5)
+ singlelimitlist.append(singlemin.Clone())
+
+ singlemax.SetLineColor(gColorList[i])
+ singlemax.SetLineWidth(2)
+ singlemax.SetLineStyle(5)
+ singlelimitlist.append(singlemax.Clone())
+
+
+ # Assign Ranges
+ miny = 99999.9
+ maxy = 0.0
+
+ for singlemc in singlemclist:
+ miny = min([singlemc.GetMinimum(),miny])
+ maxy = max([singlemc.GetMaximum(),maxy])
+ for singlelimit in singlelimitlist:
+ miny = min([singlelimit.GetMinimum(),miny])
+ maxy = max([singlelimit.GetMaximum(),maxy])
+
+ widthy = maxy - miny
+
+ # Assign Ranges to data
+ data = singlemclist[0]
+ data.GetYaxis().SetRangeUser( miny - 0.1*widthy, maxy + 0.3*widthy)
+
+ # Draw our limits
+ for i, mc in enumerate(singlemclist):
+ if i == 0: mc.Draw("E2")
+ else: mc.Draw("SAME E2")
+
+ leg = gPad.BuildLegend(0.7,0.8,1.0,1.0)
+
+
+ for limit in singlelimitlist:
+ limit.Draw("SAME HIST")
+
+ for mc in singlemclist:
+ mc.Draw("SAME E2")
+
+ startline = TLine(0.6,0.6,0.8,0.8)
+ limitline = TLine(0.6,0.6,0.8,0.8)
+ startline.SetLineStyle(7)
+ limitline.SetLineStyle(5)
+ leg.AddEntry(startline, "Start", "l")
+ leg.AddEntry(limitline, "Limits", "l")
+
+ linehists = []
+ for mc in singlemclist:
+ mcline = mc.Clone()
+ mcline.SetFillStyle(0)
+ mcline.Draw("SAME HIST")
+ linehists.append(mcline)
+
+ leg.Draw("SAME")
+ gPad.Update()
+
+
+ return [singlemclist, singlelimitlist, linehists, leg]
+
+
+
+
+
+if __name__=="__main__":
+ c1 = TCanvas("c1","c1",800,600)
+ c1.cd()
+
+
+ # Make filelist
+ allfiles = []
+
+ for i in xrange(2, len(sys.argv)):
+ print "Reading ", i, sys.argv[i]
+
+ # split by comma
+ splitname = sys.argv[i].split(",")
+
+ # Get First
+ if (os.path.isfile(splitname[0])):
+
+ # Get File
+ newfile = (TFile(splitname[0],"READ"))
+ if not newfile:
+ print "File is not a ROOT file : ", splitname[0]
+ sys.exit()
+
+ # Set Name
+ name = splitname[0].replace(".root","")
+ if len(splitname) > 1:
+ name = splitname[1]
+
+ allfiles.append([newfile, name])
+
+
+ print allfiles
+
+ # Parse Unique Keys
+ uniquekeys = []
+ for readfile in allfiles:
+ for readkey in readfile[0].GetListOfKeys():
+ if not (readkey.GetName().endswith("_data")): continue
+ if readkey.GetName() not in uniquekeys: uniquekeys.append(readkey.GetName())
+
+ print uniquekeys
+
+ # Setup First Page
+ leg = TLegend(0.1,0.1,0.9,0.9)
+ for i, readfile in enumerate(allfiles):
+ hist = TH1D(readfile[1],readfile[1],1,0,1)
+ hist.SetLineColor(gColorList[i % len(gColorList)])
+ hist.SetLineWidth(2)
+
+ leg.AddEntry(hist, readfile[1], "l")
+
+ leg.Draw()
+ gPad.Update()
+
+
+ outputfile = sys.argv[1]
+ c1.Print(outputfile + "(")
+
+ plotlist = DrawFitDialsPlot(allfiles)
+ c1.Print(outputfile)
+
+ # Loop through unique keys
+ for readkey in uniquekeys:
+
+ # Draw
+ datamclist = DrawDataMC(readkey, allfiles)
+
+ # Save
+ c1.Print(outputfile)
+
+ # Now save the legend again to close...
+ leg.Draw()
+ gPad.Update()
+ gPad.Print(outputfile + ")")
diff --git a/src/InputHandler/FitEvent.cxx b/src/InputHandler/FitEvent.cxx
index cbdf359..0d8503a 100644
--- a/src/InputHandler/FitEvent.cxx
+++ b/src/InputHandler/FitEvent.cxx
@@ -1,445 +1,448 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is pddrt of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "FitEvent.h"
#include <iostream>
#include "TObjArray.h"
FitEvent::FitEvent() {
fGenInfo = NULL;
kRemoveFSIParticles = true;
kRemoveUndefParticles = true;
AllocateParticleStack(400);
};
void FitEvent::AddGeneratorInfo(GeneratorInfoBase* gen) {
fGenInfo = gen;
gen->AllocateParticleStack(kMaxParticles);
}
void FitEvent::AllocateParticleStack(int stacksize) {
LOG(DEB) << "Allocating particle stack of size: " << stacksize << std::endl;
kMaxParticles = stacksize;
fParticleList = new FitParticle*[kMaxParticles];
fParticleMom = new double*[kMaxParticles];
fParticleState = new UInt_t[kMaxParticles];
fParticlePDG = new int[kMaxParticles];
fOrigParticleMom = new double*[kMaxParticles];
fOrigParticleState = new UInt_t[kMaxParticles];
fOrigParticlePDG = new int[kMaxParticles];
for (size_t i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
fParticleMom[i] = new double[4];
fOrigParticleMom[i] = new double[4];
}
if (fGenInfo) fGenInfo->AllocateParticleStack(kMaxParticles);
}
void FitEvent::ExpandParticleStack(int stacksize) {
DeallocateParticleStack();
AllocateParticleStack(stacksize);
}
void FitEvent::DeallocateParticleStack() {
for (size_t i = 0; i < kMaxParticles; i++) {
if (fParticleList[i]) delete fParticleList[i];
delete fParticleMom[i];
delete fOrigParticleMom[i];
}
delete fParticleMom;
delete fOrigParticleMom;
delete fParticleList;
delete fParticleState;
delete fParticlePDG;
delete fOrigParticleState;
delete fOrigParticlePDG;
if (fGenInfo) fGenInfo -> DeallocateParticleStack();
kMaxParticles = 0;
}
void FitEvent::ClearFitParticles() {
for (size_t i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
}
}
void FitEvent::FreeFitParticles() {
for (size_t i = 0; i < kMaxParticles; i++) {
FitParticle* fp = fParticleList[i];
if (fp) delete fp;
fParticleList[i] = NULL;
}
}
void FitEvent::ResetParticleList() {
for (unsigned int i = 0; i < kMaxParticles; i++) {
FitParticle* fp = fParticleList[i];
if (fp) delete fp;
fParticleList[i] = NULL;
}
}
void FitEvent::HardReset() {
for (unsigned int i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
}
}
void FitEvent::ResetEvent() {
fMode = 9999;
Mode = 9999;
fEventNo = -1;
fTotCrs = -1.0;
fTargetA = -1;
fTargetZ = -1;
fTargetH = -1;
fBound = false;
fNParticles = 0;
if (fGenInfo) fGenInfo->Reset();
for (unsigned int i = 0; i < kMaxParticles; i++) {
if (fParticleList[i]) delete fParticleList[i];
fParticleList[i] = NULL;
continue;
fParticlePDG[i] = 0;
fParticleState[i] = kUndefinedState;
fParticleMom[i][0] = 0.0;
fParticleMom[i][1] = 0.0;
fParticleMom[i][2] = 0.0;
fParticleMom[i][3] = 0.0;
fOrigParticlePDG[i] = 0;
fOrigParticleState[i] = kUndefinedState;
fOrigParticleMom[i][0] = 0.0;
fOrigParticleMom[i][1] = 0.0;
fOrigParticleMom[i][2] = 0.0;
fOrigParticleMom[i][3] = 0.0;
}
}
void FitEvent::OrderStack() {
// Copy current stack
int npart = fNParticles;
for (int i = 0; i < npart; i++) {
fOrigParticlePDG[i] = fParticlePDG[i];
fOrigParticleState[i] = fParticleState[i];
fOrigParticleMom[i][0] = fParticleMom[i][0];
fOrigParticleMom[i][1] = fParticleMom[i][1];
fOrigParticleMom[i][2] = fParticleMom[i][2];
fOrigParticleMom[i][3] = fParticleMom[i][3];
}
// Now run loops for each particle
fNParticles = 0;
int stateorder[6] = {kInitialState, kFinalState, kFSIState,
kNuclearInitial, kNuclearRemnant, kUndefinedState
};
for (int s = 0; s < 6; s++) {
for (int i = 0; i < npart; i++) {
if ((UInt_t)fOrigParticleState[i] != (UInt_t)stateorder[s]) continue;
fParticlePDG[fNParticles] = fOrigParticlePDG[i];
fParticleState[fNParticles] = fOrigParticleState[i];
fParticleMom[fNParticles][0] = fOrigParticleMom[i][0];
fParticleMom[fNParticles][1] = fOrigParticleMom[i][1];
fParticleMom[fNParticles][2] = fOrigParticleMom[i][2];
fParticleMom[fNParticles][3] = fOrigParticleMom[i][3];
fNParticles++;
}
}
if (LOG_LEVEL(DEB)) {
LOG(DEB) << "Ordered stack" << std::endl;
for (int i = 0; i < fNParticles; i++) {
LOG(DEB) << "Particle " << i << ". " << fParticlePDG[i] << " "
<< fParticleMom[i][0] << " " << fParticleMom[i][1] << " "
<< fParticleMom[i][2] << " " << fParticleMom[i][3] << " "
<< fParticleState[i] << std::endl;
}
}
if (fNParticles != npart) {
ERR(FTL) << "Dropped some particles when ordering the stack!" << std::endl;
}
return;
}
void FitEvent::Print() {
if (LOG_LEVEL(EVT)) {
LOG(EVT) << "EVTEvent print" << std::endl;
LOG(EVT) << "Mode: " << fMode << std::endl;
LOG(EVT) << "Particles: " << fNParticles << std::endl;
LOG(EVT) << " -> Particle Stack " << std::endl;
for (int i = 0; i < fNParticles; i++) {
LOG(EVT) << " -> -> " << i << ". " << fParticlePDG[i] << " "
<< fParticleState[i] << " "
<< " Mom(" << fParticleMom[i][0] << ", " << fParticleMom[i][1]
<< ", " << fParticleMom[i][2] << ", " << fParticleMom[i][3] << ")."
<< std::endl;
}
}
return;
}
/* Read/Write own event class */
void FitEvent::SetBranchAddress(TChain* tn) {
tn->SetBranchAddress("Mode", &fMode);
tn->SetBranchAddress("Mode", &Mode);
tn->SetBranchAddress("EventNo", &fEventNo);
tn->SetBranchAddress("TotCrs", &fTotCrs);
tn->SetBranchAddress("TargetA", &fTargetA);
tn->SetBranchAddress("TargetH", &fTargetH);
tn->SetBranchAddress("Bound", &fBound);
tn->SetBranchAddress("RWWeight", &SavedRWWeight);
tn->SetBranchAddress("InputWeight", &InputWeight);
- tn->SetBranchAddress("NParticles", &fNParticles);
- tn->SetBranchAddress("ParticleState", fParticleState);
- tn->SetBranchAddress("ParticlePDG", fParticlePDG);
- tn->SetBranchAddress("ParticleMom", fParticleMom);
+
+
+ // This has to be setup by the handler now :(
+ // tn->SetBranchAddress("NParticles", &fNParticles);
+ // tn->SetBranchAddress("ParticleState", &fParticleState);
+ // tn->SetBranchAddress("ParticlePDG", &fParticlePDG);
+ // tn->SetBranchAddress("ParticleMom", &fParticleMom);
}
void FitEvent::AddBranchesToTree(TTree* tn) {
tn->Branch("Mode", &Mode, "Mode/I");
tn->Branch("EventNo", &fEventNo, "EventNo/i");
tn->Branch("TotCrs", &fTotCrs, "TotCrs/D");
tn->Branch("TargetA", &fTargetA, "TargetA/I");
tn->Branch("TargetH", &fTargetH, "TargetH/I");
tn->Branch("Bound", &fBound, "Bound/O");
tn->Branch("RWWeight", &RWWeight, "RWWeight/D");
tn->Branch("InputWeight", &InputWeight, "InputWeight/D");
tn->Branch("NParticles", &fNParticles, "NParticles/I");
tn->Branch("ParticleState", fOrigParticleState, "ParticleState[NParticles]/i");
tn->Branch("ParticlePDG", fOrigParticlePDG, "ParticlePDG[NParticles]/I");
tn->Branch("ParticleMom", fOrigParticleMom, "ParticleMom[NParticles][4]/D");
}
// ------- EVENT ACCESS FUNCTION --------- //
TLorentzVector FitEvent::GetParticleP4 (int index) const {
if (index == -1 or index >= fNParticles) return TLorentzVector();
return TLorentzVector( fParticleMom[index][0],
fParticleMom[index][1],
fParticleMom[index][2],
fParticleMom[index][3] );
}
TVector3 FitEvent::GetParticleP3 (int index) const {
if (index == -1 or index >= fNParticles) return TVector3();
return TVector3( fParticleMom[index][0],
fParticleMom[index][1],
fParticleMom[index][2] );
}
double FitEvent::GetParticleMom (int index) const {
if (index == -1 or index >= fNParticles) return 0.0;
return sqrt(fParticleMom[index][0] * fParticleMom[index][0] +
fParticleMom[index][1] * fParticleMom[index][1] +
fParticleMom[index][2] * fParticleMom[index][2]);
}
double FitEvent::GetParticleMom2 (int index) const {
if (index == -1 or index >= fNParticles) return 0.0;
return fabs((fParticleMom[index][0] * fParticleMom[index][0] +
fParticleMom[index][1] * fParticleMom[index][1] +
fParticleMom[index][2] * fParticleMom[index][2]));
}
double FitEvent::GetParticleE (int index) const {
if (index == -1 or index >= fNParticles) return 0.0;
return fParticleMom[index][3];
}
int FitEvent::GetParticleState (int index) const {
if (index == -1 or index >= fNParticles) return kUndefinedState;
return (fParticleState[index]);
}
int FitEvent::GetParticlePDG (int index) const {
if (index == -1 or index >= fNParticles) return 0;
return (fParticlePDG[index]);
}
FitParticle* FitEvent::GetParticle (int const i) {
// Check Valid Index
if (i == -1){
return NULL;
}
// Check Valid
if (i > fNParticles) {
ERR(FTL) << "Requesting particle beyond stack!" << std::endl
<< "i = " << i << " N = " << fNParticles << std::endl
<< "Mode = " << fMode << std::endl;
throw;
}
if (!fParticleList[i]) {
fParticleList[i] = new FitParticle(fParticleMom[i][0], fParticleMom[i][1],
fParticleMom[i][2], fParticleMom[i][3],
fParticlePDG[i], fParticleState[i]);
} else {
fParticleList[i]->SetValues(fParticleMom[i][0], fParticleMom[i][1],
fParticleMom[i][2], fParticleMom[i][3],
fParticlePDG[i], fParticleState[i]);
}
return fParticleList[i];
}
bool FitEvent::HasParticle(int const pdg, int const state) const {
bool found = false;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 && fParticleState[i] != (uint)state) continue;
if (fParticlePDG[i] == pdg) found = true;
}
return found;
}
int FitEvent::NumParticle(int const pdg, int const state) const {
int nfound = 0;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state) continue;
if (pdg == 0 or fParticlePDG[i] == pdg) nfound += 1;
}
return nfound;
}
std::vector<int> FitEvent::GetAllParticleIndices (int const pdg, int const state) const {
std::vector<int> indexlist;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state) continue;
if (pdg == 0 or fParticlePDG[i] == pdg) {
indexlist.push_back(i);
}
}
return indexlist;
}
std::vector<FitParticle*> FitEvent::GetAllParticle(int const pdg, int const state) {
std::vector<int> indexlist = GetAllParticleIndices(pdg, state);
std::vector<FitParticle*> plist;
for (std::vector<int>::iterator iter = indexlist.begin();
iter != indexlist.end(); iter++) {
plist.push_back( GetParticle((*iter)) );
}
return plist;
}
int FitEvent::GetHMParticleIndex (int const pdg, int const state) const {
double maxmom2 = -9999999.9;
int maxind = -1;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state) continue;
if (pdg == 0 or fParticlePDG[i] == pdg) {
double newmom2 = GetParticleMom2(i);
if (newmom2 > maxmom2) {
maxind = i;
maxmom2 = newmom2;
}
}
}
return maxind;
}
int FitEvent::GetBeamNeutrinoIndex (void) const {
for (int i = 0; i < fNParticles; i++){
if (fParticleState[i] != kInitialState) continue;
int pdg = abs(fParticlePDG[i]);
if (pdg == 12 or pdg == 14 or pdg == 16){
return i;
}
}
return 0;
}
int FitEvent::GetBeamElectronIndex (void) const {
return GetHMISParticleIndex( 11 );
}
int FitEvent::GetBeamPionIndex (void) const {
return GetHMISParticleIndex( PhysConst::pdg_pions );
}
int FitEvent::NumFSMesons() {
int nMesons = 0;
for (int i = 0; i < fNParticles; i++) {
if (fParticleState[i] != kFinalState) continue;
if (abs(fParticlePDG[i]) >= 111 && abs(fParticlePDG[i]) <= 557)
nMesons += 1;
}
return nMesons;
}
int FitEvent::NumFSLeptons(void) const{
int nLeptons = 0;
for (int i = 0; i < fNParticles; i++) {
if (fParticleState[i] != kFinalState) continue;
if (abs(fParticlePDG[i]) == 11 || abs(fParticlePDG[i]) == 13 ||
abs(fParticlePDG[i]) == 15)
nLeptons += 1;
}
return nLeptons;
}
diff --git a/src/InputHandler/FitEventInputHandler.cxx b/src/InputHandler/FitEventInputHandler.cxx
index 7de683a..969c245 100644
--- a/src/InputHandler/FitEventInputHandler.cxx
+++ b/src/InputHandler/FitEventInputHandler.cxx
@@ -1,100 +1,130 @@
#include "FitEventInputHandler.h"
FitEventInputHandler::FitEventInputHandler(std::string const& handle, std::string const& rawinputs) {
LOG(SAM) << "Creating FitEventInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
fFitEventTree = new TChain("nuisance_events");
fCacheSize = FitPar::Config().GetParI("CacheSize");
std::vector<std::string> inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile* inp_file = new TFile(inputs[inp_it].c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
ERR(FTL) << "FitEvent File IsZombie() at " << inputs[inp_it] << std::endl;
throw;
}
// Get Flux/Event hist
TH1D* fluxhist = (TH1D*)inp_file->Get("nuisance_fluxhist");
TH1D* eventhist = (TH1D*)inp_file->Get("nuisance_eventhist");
if (!fluxhist or !eventhist) {
ERR(FTL) << "FitEvent FILE doesn't contain flux/xsec info" << std::endl;
throw;
}
// Get N Events
TTree* eventtree = (TTree*)inp_file->Get("nuisance_events");
if (!eventtree) {
ERR(FTL) << "nuisance_events not located in GENIE file! " << inputs[inp_it] << std::endl;
throw;
}
int nevents = eventtree->GetEntries();
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add to TChain
fFitEventTree->Add( inputs[inp_it].c_str() );
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kINPUTFITEVENT;
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->HardReset();
fNUISANCEEvent->SetBranchAddress(fFitEventTree);
-
+
+ fFitEventTree->SetBranchAddress("NParticles", &fReadNParticles);
+ fFitEventTree->SetBranchAddress("ParticleState", &fReadParticleState);
+ fFitEventTree->SetBranchAddress("ParticlePDG", &fReadParticlePDG);
+ fFitEventTree->SetBranchAddress("ParticleMom", &fReadParticleMom);
+
+
+ fFitEventTree->Show(0);
+ fNUISANCEEvent = GetNuisanceEvent(0);
+ std::cout << "NParticles = " << fNUISANCEEvent->Npart() << std::endl;
+ std::cout << "Event Info " << fNUISANCEEvent->PartInfo(0)->fPID << std::endl;
+
}
FitEventInputHandler::~FitEventInputHandler(){
if (fFitEventTree) delete fFitEventTree;
}
void FitEventInputHandler::CreateCache() {
- fFitEventTree->SetCacheEntryRange(0, fNEvents);
- fFitEventTree->AddBranchToCache("*", 1);
- fFitEventTree->SetCacheSize(fCacheSize);
+ // fFitEventTree->SetCacheEntryRange(0, fNEvents);
+ // fFitEventTree->AddBranchToCache("*", 1);
+ // fFitEventTree->SetCacheSize(fCacheSize);
}
void FitEventInputHandler::RemoveCache() {
- fFitEventTree->SetCacheEntryRange(0, fNEvents);
- fFitEventTree->AddBranchToCache("*", 0);
- fFitEventTree->SetCacheSize(0);
+ //fFitEventTree->SetCacheEntryRange(0, fNEvents);
+ // fFitEventTree->AddBranchToCache("*", 0);
+ // fFitEventTree->SetCacheSize(0);
}
FitEvent* FitEventInputHandler::GetNuisanceEvent(const UInt_t entry, const bool lightweight) {
// Return NULL if out of bounds
if (entry >= (UInt_t)fNEvents) return NULL;
// Reset all variables before tree read
fNUISANCEEvent->ResetEvent();
// Read NUISANCE Tree
fFitEventTree->GetEntry(entry);
+ // Fill Stack
+ fNUISANCEEvent->fNParticles = 0;
+ for (size_t i = 0; i < fReadNParticles; i++){
+ size_t curpart = fNUISANCEEvent->fNParticles;
+ fNUISANCEEvent->fParticleState[curpart] = fReadParticleState[i];
+
+ // Mom
+ fNUISANCEEvent->fParticleMom[curpart][0] = fReadParticleMom[i][0];
+ fNUISANCEEvent->fParticleMom[curpart][1] = fReadParticleMom[i][1];
+ fNUISANCEEvent->fParticleMom[curpart][2] = fReadParticleMom[i][2];
+ fNUISANCEEvent->fParticleMom[curpart][3] = fReadParticleMom[i][3];
+
+ // PDG
+ fNUISANCEEvent->fParticlePDG[curpart] = fReadParticlePDG[i];
+
+ // Add to N particle count
+ fNUISANCEEvent->fNParticles++;
+ }
+
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
return fNUISANCEEvent;
}
double FitEventInputHandler::GetInputWeight(int entry) {
double w = InputHandlerBase::GetInputWeight(entry);
return w * fNUISANCEEvent->SavedRWWeight;
}
void FitEventInputHandler::Print() {
}
diff --git a/src/InputHandler/FitEventInputHandler.h b/src/InputHandler/FitEventInputHandler.h
index 45eb979..5bdc84c 100644
--- a/src/InputHandler/FitEventInputHandler.h
+++ b/src/InputHandler/FitEventInputHandler.h
@@ -1,37 +1,43 @@
#ifndef FITEVENT_INPUTHANDLER_H
#define FITEVENT_INPUTHANDLER_H
/*!
* \addtogroup InputHandler
* @{
*/
#include "InputHandler.h"
#include "FitEvent.h"
#include "PlotUtils.h"
/// Class to read in NUISANCE FitEvents that have been saved to tree
class FitEventInputHandler : public InputHandlerBase {
public:
/// Standard constructor given name and inputs
FitEventInputHandler(std::string const& handle, std::string const& rawinputs);
virtual ~FitEventInputHandler();
/// Create a TTree Cache to speed up file read
void CreateCache();
/// Remove TTree Cache to save memory
void RemoveCache();
/// Returns NUISANCE FitEvent from the TTree. If lightweight does nothing.
FitEvent* GetNuisanceEvent(const UInt_t entry, const bool lightweight=false);
/// Alongside InputWeight also returns any saved RWWeights
double GetInputWeight(int entry);
/// Print out event information
void Print();
TChain* fFitEventTree; ///< TTree from FitEvent file.
+
+ int fReadNParticles;
+ double fReadParticleMom[400][4];
+ UInt_t fReadParticleState[400];
+ int fReadParticlePDG[400];
+
};
/*! @} */
#endif
diff --git a/src/InputHandler/InputHandler.cxx b/src/InputHandler/InputHandler.cxx
index b11a028..c21c968 100644
--- a/src/InputHandler/InputHandler.cxx
+++ b/src/InputHandler/InputHandler.cxx
@@ -1,263 +1,263 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "InputHandler.h"
InputHandlerBase::InputHandlerBase() {
fName = "";
fFluxHist = NULL;
fEventHist = NULL;
fNEvents = 0;
fNUISANCEEvent = NULL;
fBaseEvent = NULL;
kRemoveUndefParticles = FitPar::Config().GetParB("RemoveUndefParticles");
kRemoveFSIParticles = FitPar::Config().GetParB("RemoveFSIParticles");
kRemoveNuclearParticles = FitPar::Config().GetParB("RemoveNuclearParticles");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
fTTreePerformance = NULL;
};
InputHandlerBase::~InputHandlerBase() {
if (fFluxHist) delete fFluxHist;
if (fEventHist) delete fEventHist;
- if (fXSecHist) delete fXSecHist;
- if (fNUISANCEEvent) delete fNUISANCEEvent;
+ // if (fXSecHist) delete fXSecHist;
+ // if (fNUISANCEEvent) delete fNUISANCEEvent;
jointfluxinputs.clear();
jointeventinputs.clear();
jointindexlow.clear();
jointindexhigh.clear();
jointindexallowed.clear();
jointindexscale.clear();
- if (fTTreePerformance) {
- fTTreePerformance->SaveAs(("ttreeperfstats_" + fName + ".root").c_str());
- }
+ // if (fTTreePerformance) {
+ // fTTreePerformance->SaveAs(("ttreeperfstats_" + fName + ".root").c_str());
+ // }
}
void InputHandlerBase::Print() {
};
TH1D* InputHandlerBase::GetXSecHistogram(void) {
fXSecHist = (TH1D*)fFluxHist->Clone();
fXSecHist->Divide(fEventHist);
return fXSecHist;
};
double InputHandlerBase::PredictedEventRate(double low, double high,
std::string intOpt) {
int minBin = fFluxHist->GetXaxis()->FindBin(low);
int maxBin = fFluxHist->GetXaxis()->FindBin(high);
return fEventHist->Integral(minBin, maxBin + 1, intOpt.c_str());
};
double InputHandlerBase::TotalIntegratedFlux(double low, double high,
std::string intOpt) {
Int_t minBin = fFluxHist->GetXaxis()->FindFixBin(low);
Int_t maxBin = fFluxHist->GetXaxis()->FindFixBin(high);
if ((fFluxHist->IsBinOverflow(minBin) && (low != -9999.9))) {
minBin = 1;
}
if ((fFluxHist->IsBinOverflow(maxBin) && (high != -9999.9))) {
maxBin = fFluxHist->GetXaxis()->GetNbins() + 1;
}
// If we are within a single bin
if (minBin == maxBin) {
// Get the contained fraction of the single bin's width
return ((high - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) *
fFluxHist->Integral(minBin, minBin, intOpt.c_str());
}
double lowBinUpEdge = fFluxHist->GetXaxis()->GetBinUpEdge(minBin);
double highBinLowEdge = fFluxHist->GetXaxis()->GetBinLowEdge(maxBin);
double lowBinfracIntegral =
((lowBinUpEdge - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) *
fFluxHist->Integral(minBin, minBin, intOpt.c_str());
double highBinfracIntegral =
((high - highBinLowEdge) / fFluxHist->GetXaxis()->GetBinWidth(maxBin)) *
fFluxHist->Integral(maxBin, maxBin, intOpt.c_str());
// If they are neighbouring bins
if ((minBin + 1) == maxBin) {
std::cout << "Get lowfrac + highfrac" << std::endl;
// Get the contained fraction of the two bin's width
return lowBinfracIntegral + highBinfracIntegral;
}
// If there are filled bins between them
return lowBinfracIntegral + highBinfracIntegral +
fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str());
// return fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str());
}
std::vector<TH1*> InputHandlerBase::GetFluxList(void) {
return std::vector<TH1*>(1, fFluxHist);
};
std::vector<TH1*> InputHandlerBase::GetEventList(void) {
return std::vector<TH1*>(1, fEventHist);
};
std::vector<TH1*> InputHandlerBase::GetXSecList(void) {
return std::vector<TH1*>(1, GetXSecHistogram());
};
FitEvent* InputHandlerBase::FirstNuisanceEvent() {
fCurrentIndex = 0;
return GetNuisanceEvent(fCurrentIndex);
};
FitEvent* InputHandlerBase::NextNuisanceEvent() {
fCurrentIndex++;
return GetNuisanceEvent(fCurrentIndex);
};
BaseFitEvt* InputHandlerBase::FirstBaseEvent() {
fCurrentIndex = 0;
return GetBaseEvent(fCurrentIndex);
};
BaseFitEvt* InputHandlerBase::NextBaseEvent() {
fCurrentIndex++;
if (jointinput and fMaxEvents != -1) {
while ( fCurrentIndex < jointindexlow[jointindexswitch] ||
fCurrentIndex >= jointindexhigh[jointindexswitch] ) {
jointindexswitch++;
// Loop Around
if (jointindexswitch == jointindexlow.size()) {
jointindexswitch = 0;
}
}
if (fCurrentIndex > jointindexlow[jointindexswitch] + jointindexallowed[jointindexswitch]) {
fCurrentIndex = jointindexlow[jointindexswitch];
}
}
return GetBaseEvent(fCurrentIndex);
};
void InputHandlerBase::RegisterJointInput(std::string input, int n, TH1D* f, TH1D* e) {
if (jointfluxinputs.size() == 0) {
jointindexswitch = 0;
fNEvents = 0;
}
// Push into individual input vectors
jointfluxinputs.push_back( (TH1D*) f->Clone() );
jointeventinputs.push_back( (TH1D*) e->Clone() );
jointindexlow.push_back(fNEvents);
jointindexhigh.push_back(fNEvents + n);
fNEvents += n;
// Add to the total flux/event hist
if (!fFluxHist) fFluxHist = (TH1D*) f->Clone();
else fFluxHist->Add(f);
if (!fEventHist) fEventHist = (TH1D*) e->Clone();
else fEventHist->Add(e);
}
void InputHandlerBase::SetupJointInputs() {
if (jointeventinputs.size() <= 1) {
jointinput = false;
} else if (jointeventinputs.size() > 1) {
jointinput = true;
jointindexswitch = 0;
}
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
if (fMaxEvents != -1 and jointeventinputs.size() > 1){
THROW("Can only handle joint inputs when config MAXEVENTS = -1!");
}
for (size_t i = 0; i < jointeventinputs.size(); i++) {
TH1D* eventhist = (TH1D*) jointeventinputs.at(i)->Clone();
double scale = double(fNEvents) / fEventHist->Integral("width");
scale *= eventhist->Integral("width");
scale /= double(jointindexhigh[i] - jointindexlow[i]);
jointindexscale .push_back(scale);
}
fEventHist->SetNameTitle((fName + "_EVT").c_str(), (fName + "_EVT").c_str());
fFluxHist->SetNameTitle((fName + "_FLUX").c_str(), (fName + "_FLUX").c_str());
// Setup Max Events
if (fMaxEvents > 1 && fMaxEvents < fNEvents) {
if (LOG_LEVEL(SAM)) {
std::cout << "\t\t|-> Read Max Entries : " << fMaxEvents << std::endl;
}
fNEvents = fMaxEvents;
}
// Print out Status
if (LOG_LEVEL(SAM)) {
std::cout << "\t\t|-> Total Entries : " << fNEvents << std::endl
<< "\t\t|-> Event Integral : " << fEventHist->Integral("width") * 1.E-38 << " events/nucleon" << std::endl
<< "\t\t|-> Flux Integral : " << fFluxHist->Integral("width") << " /cm2" << std::endl
<< "\t\t|-> Event/Flux : "
<< fEventHist->Integral("width") * 1.E-38 / fFluxHist->Integral("width") << " cm2/nucleon" << std::endl;
}
}
BaseFitEvt* InputHandlerBase::GetBaseEvent(const UInt_t entry) {
return static_cast<BaseFitEvt*>(GetNuisanceEvent(entry, true));
}
double InputHandlerBase::GetInputWeight(int entry) {
if (!jointinput) return 1.0;
// Find Switch Scale
while ( entry < jointindexlow[jointindexswitch] ||
entry >= jointindexhigh[jointindexswitch] ) {
jointindexswitch++;
// Loop Around
if (jointindexswitch >= jointindexlow.size()) {
jointindexswitch = 0;
}
}
return jointindexscale[jointindexswitch];
};
diff --git a/src/InputHandler/NEUTInputHandler.cxx b/src/InputHandler/NEUTInputHandler.cxx
index fb46488..8883c8e 100644
--- a/src/InputHandler/NEUTInputHandler.cxx
+++ b/src/InputHandler/NEUTInputHandler.cxx
@@ -1,436 +1,436 @@
#ifdef __NEUT_ENABLED__
#include "NEUTInputHandler.h"
NEUTGeneratorInfo::~NEUTGeneratorInfo() {
DeallocateParticleStack();
}
void NEUTGeneratorInfo::AddBranchesToTree(TTree * tn) {
tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I");
tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode, "NEUTParticleStatusCode[NEUTParticleN]/I");
tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode, "NEUTParticleAliveCode[NEUTParticleN]/I");
}
void NEUTGeneratorInfo::SetBranchesFromTree(TTree* tn) {
tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN );
tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode );
tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode );
}
void NEUTGeneratorInfo::AllocateParticleStack(int stacksize) {
fNEUTParticleN = 0;
fNEUTParticleStatusCode = new int[stacksize];
fNEUTParticleStatusCode = new int[stacksize];
}
void NEUTGeneratorInfo::DeallocateParticleStack() {
delete fNEUTParticleStatusCode;
delete fNEUTParticleAliveCode;
}
void NEUTGeneratorInfo::FillGeneratorInfo(NeutVect* nevent) {
Reset();
for (int i = 0; i < nevent->Npart(); i++) {
fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus;
fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive;
fNEUTParticleN++;
}
}
void NEUTGeneratorInfo::Reset() {
for (int i = 0; i < fNEUTParticleN; i++) {
fNEUTParticleStatusCode[i] = -1;
fNEUTParticleAliveCode[i] = 9;
}
fNEUTParticleN = 0;
}
NEUTInputHandler::NEUTInputHandler(std::string const& handle, std::string const& rawinputs) {
LOG(SAM) << "Creating NEUTInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
// Setup the TChain
fNEUTTree = new TChain("neuttree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraNEUT");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector<std::string> inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile* inp_file = new TFile(inputs[inp_it].c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
THROW( "NEUT File IsZombie() at : '" << inputs[inp_it] << "'" << std::endl
<< "Check that your file paths are correct and the file exists!" << std::endl
<< "$ ls -lh " << inputs[inp_it] );
}
// Get Flux/Event hist
TH1D* fluxhist = (TH1D*)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "flux")).c_str());
TH1D* eventhist = (TH1D*)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "evt")).c_str());
if (!fluxhist or !eventhist) {
ERROR(FTL, "Input File Contents: " << inputs[inp_it] );
inp_file->ls();
THROW( "NEUT FILE doesn't contain flux/xsec info. You may have to regenerate your MC!" );
}
// Get N Events
TTree* neuttree = (TTree*)inp_file->Get("neuttree");
if (!neuttree) {
ERROR(FTL, "neuttree not located in NEUT file: " << inputs[inp_it]);
THROW("Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = neuttree->GetEntries();
if (nevents <= 0) {
THROW("Trying to a TTree with " << nevents << " to TChain from : " << inputs[inp_it]);
}
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add To TChain
fNEUTTree->AddFile( inputs[inp_it].c_str() );
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kNEUT;
fNeutVect = NULL;
fNEUTTree->SetBranchAddress("vectorbranch", &fNeutVect);
fNEUTTree->GetEntry(0);
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetNeutVect(fNeutVect);
if (fSaveExtra) {
fNeutInfo = new NEUTGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fNeutInfo);
}
fNUISANCEEvent->HardReset();
};
NEUTInputHandler::~NEUTInputHandler() {
- if (fNEUTTree) delete fNEUTTree;
- if (fNeutVect) delete fNeutVect;
- if (fNeutInfo) delete fNeutInfo;
+ // if (fNEUTTree) delete fNEUTTree;
+ // if (fNeutVect) delete fNeutVect;
+ // if (fNeutInfo) delete fNeutInfo;
};
void NEUTInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fNEUTTree->SetCacheEntryRange(0, fNEvents);
fNEUTTree->AddBranchToCache("vectorbranch", 1);
fNEUTTree->SetCacheSize(fCacheSize);
}
}
void NEUTInputHandler::RemoveCache() {
// fNEUTTree->SetCacheEntryRange(0, fNEvents);
fNEUTTree->AddBranchToCache("vectorbranch", 0);
fNEUTTree->SetCacheSize(0);
}
FitEvent* NEUTInputHandler::GetNuisanceEvent(const UInt_t entry, const bool lightweight) {
// Catch too large entries
if (entry >= (UInt_t)fNEvents) return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fNEUTTree->GetEntry(entry);
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
// Return event pointer
return fNUISANCEEvent;
}
int NEUTInputHandler::GetNeutParticleStatus(NeutPart * part) {
// State
int state = kUndefinedState;
// fStatus == -1 means initial state
if (part->fIsAlive == false && part->fStatus == -1) {
state = kInitialState;
// NEUT has a bit of a strange convention for fIsAlive and fStatus
// combinations
// for NC and neutrino particle isAlive true/false and status 2 means
// final state particle
// for other particles in NC status 2 means it's an FSI particle
// for CC it means it was an FSI particle
} else if (part->fStatus == 2) {
// NC case is a little strange... The outgoing neutrino might be alive or
// not alive. Remaining particles with status 2 are FSI particles that
// reinteracted
if (abs(fNeutVect->Mode) > 30 &&
(abs(part->fPID) == 14 || abs(part->fPID) == 12)) {
state = kFinalState;
// The usual CC case
} else if (part->fIsAlive == true) {
state = kFSIState;
}
} else if (part->fIsAlive == true && part->fStatus == 2 &&
(abs(part->fPID) == 14 || abs(part->fPID) == 12)) {
state = kFinalState;
} else if (part->fIsAlive == true && part->fStatus == 0) {
state = kFinalState;
} else if (part->fIsAlive == true) {
ERR(WRN) << "Undefined NEUT state "
<< " Alive: " << part->fIsAlive << " Status: " << part->fStatus
<< " PDG: " << part->fPID << std::endl;
throw;
}
return state;
}
void NEUTInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Fill Globals
fNUISANCEEvent->fMode = fNeutVect->Mode;
fNUISANCEEvent->Mode = fNeutVect->Mode;
fNUISANCEEvent->fEventNo = fNeutVect->EventNo;
fNUISANCEEvent->fTargetA = fNeutVect->TargetA;
fNUISANCEEvent->fTargetZ = fNeutVect->TargetZ;
fNUISANCEEvent->fTargetH = fNeutVect->TargetH;
fNUISANCEEvent->fBound = bool(fNeutVect->Ibound);
if (fNUISANCEEvent->fBound) {
fNUISANCEEvent->fTargetPDG = TargetUtils::GetTargetPDGFromZA(fNUISANCEEvent->fTargetZ,
fNUISANCEEvent->fTargetA);
} else {
fNUISANCEEvent->fTargetPDG = 1000010010;
}
// Check Particle Stack
UInt_t npart = fNeutVect->Npart();
UInt_t kmax = fNUISANCEEvent->kMaxParticles;
if (npart > kmax) {
ERR(FTL) << "NEUT has too many particles. Expanding stack." << std::endl;
fNUISANCEEvent->ExpandParticleStack(npart);
throw;
}
// Fill Particle Stack
for (size_t i = 0; i < npart; i++) {
// Get Current Count
int curpart = fNUISANCEEvent->fNParticles;
// Get NEUT Particle
NeutPart* part = fNeutVect->PartInfo(i);
// State
int state = GetNeutParticleStatus(part);
// Remove Undefined
if (kRemoveUndefParticles &&
state == kUndefinedState) continue;
// Remove FSI
if (kRemoveFSIParticles &&
state == kFSIState) continue;
// Remove Nuclear
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant)) continue;
// State
fNUISANCEEvent->fParticleState[curpart] = state;
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = part->fP.X();
fNUISANCEEvent->fParticleMom[curpart][1] = part->fP.Y();
fNUISANCEEvent->fParticleMom[curpart][2] = part->fP.Z();
fNUISANCEEvent->fParticleMom[curpart][3] = part->fP.T();
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = part->fPID;
// Add up particle count
fNUISANCEEvent->fNParticles++;
}
// Save Extra Generator Info
if (fSaveExtra) {
fNeutInfo->FillGeneratorInfo(fNeutVect);
}
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent-> OrderStack();
return;
}
void NEUTUtils::FillNeutCommons(NeutVect* nvect) {
// WARNING: This has only been implemented for a neuttree and not GENIE
// This should be kept in sync with T2KNIWGUtils::GetNIWGEvent(TTree)
//NEUT version info. Can't get it to compile properly with this yet
//neutversion_.corev = nvect->COREVer;
//neutversion_.nucev = nvect->NUCEVer;
//neutversion_.nuccv = nvect->NUCCVer;
// Documentation: See nework.h
nework_.modene = nvect->Mode;
nework_.numne = nvect->Npart();
nemdls_.mdlqeaf = nvect->QEVForm;
nemdls_.mdlqe = nvect->QEModel;
nemdls_.mdlspi = nvect->SPIModel;
nemdls_.mdldis = nvect->DISModel;
nemdls_.mdlcoh = nvect->COHModel;
neutcoh_.necohepi = nvect->COHModel;
nemdls_.xmaqe = nvect->QEMA;
nemdls_.xmvqe = nvect->QEMV;
nemdls_.kapp = nvect->KAPPA;
//nemdls_.sccfv = SCCFVdef;
//nemdls_.sccfa = SCCFAdef;
//nemdls_.fpqe = FPQEdef;
nemdls_.xmaspi = nvect->SPIMA;
nemdls_.xmvspi = nvect->SPIMV;
nemdls_.xmares = nvect->RESMA;
nemdls_.xmvres = nvect->RESMV;
neut1pi_.xmanffres = nvect->SPIMA;
neut1pi_.xmvnffres = nvect->SPIMV;
neut1pi_.xmarsres = nvect->RESMA;
neut1pi_.xmvrsres = nvect->RESMV;
neut1pi_.neiff = nvect->SPIForm;
neut1pi_.nenrtype = nvect->SPINRType;
neut1pi_.rneca5i = nvect->SPICA5I;
neut1pi_.rnebgscl = nvect->SPIBGScale;
nemdls_.xmacoh = nvect->COHMA;
nemdls_.rad0nu = nvect->COHR0;
//nemdls_.fa1coh = nvect->COHA1err;
//nemdls_.fb1coh = nvect->COHb1err;
//neutdis_.nepdf = NEPDFdef;
//neutdis_.nebodek = NEBODEKdef;
neutcard_.nefrmflg = nvect->FrmFlg;
neutcard_.nepauflg = nvect->PauFlg;
neutcard_.nenefo16 = nvect->NefO16;
neutcard_.nemodflg = nvect->ModFlg;
//neutcard_.nenefmodl = 1;
//neutcard_.nenefmodh = 1;
//neutcard_.nenefkinh = 1;
//neutpiabs_.neabspiemit = 1;
nenupr_.iformlen = nvect->FormLen;
neutpiless_.ipilessdcy = nvect->IPilessDcy;
neutpiless_.rpilessdcy = nvect->RPilessDcy;
neutpiless_.ipilessdcy = nvect->IPilessDcy;
neutpiless_.rpilessdcy = nvect->RPilessDcy;
neffpr_.fefqe = nvect->NuceffFactorPIQE;
neffpr_.fefqeh = nvect->NuceffFactorPIQEH;
neffpr_.fefinel = nvect->NuceffFactorPIInel;
neffpr_.fefabs = nvect->NuceffFactorPIAbs;
neffpr_.fefcx = nvect->NuceffFactorPICX;
neffpr_.fefcxh = nvect->NuceffFactorPICXH;
neffpr_.fefcoh = nvect->NuceffFactorPICoh;
neffpr_.fefqehf = nvect->NuceffFactorPIQEHKin;
neffpr_.fefcxhf = nvect->NuceffFactorPICXKin;
neffpr_.fefcohf = nvect->NuceffFactorPIQELKin;
for ( int i = 0; i < nework_.numne; i++ ) {
nework_.ipne[i] = nvect->PartInfo(i)->fPID;
nework_.pne[i][0] = (float)nvect->PartInfo(i)->fP.X() / 1000; // VC(NE)WORK in M(G)eV
nework_.pne[i][1] = (float)nvect->PartInfo(i)->fP.Y() / 1000; // VC(NE)WORK in M(G)eV
nework_.pne[i][2] = (float)nvect->PartInfo(i)->fP.Z() / 1000; // VC(NE)WORK in M(G)eV
}
// fsihist.h
// neutroot fills a dummy object for events with no FSI to prevent memory leak when
// reading the TTree, so check for it here
if ( (int)nvect->NfsiVert() == 1 ) { // An event with FSI must have at least two vertices
// if (nvect->NfsiPart()!=1 || nvect->Fsiprob!=-1)
// ERR(WRN) << "T2KNeutUtils::fill_neut_commons(TTree) NfsiPart!=1 or Fsiprob!=-1 when NfsiVert==1" << std::endl;
fsihist_.nvert = 0;
fsihist_.nvcvert = 0;
fsihist_.fsiprob = 1;
}
else { // Real FSI event
fsihist_.nvert = (int)nvect->NfsiVert();
for (int ivert = 0; ivert < fsihist_.nvert; ivert++) {
fsihist_.iflgvert[ivert] = nvect->FsiVertInfo(ivert)->fVertID;
fsihist_.posvert[ivert][0] = (float)nvect->FsiVertInfo(ivert)->fPos.X();
fsihist_.posvert[ivert][1] = (float)nvect->FsiVertInfo(ivert)->fPos.Y();
fsihist_.posvert[ivert][2] = (float)nvect->FsiVertInfo(ivert)->fPos.Z();
}
fsihist_.nvcvert = nvect->NfsiPart();
for (int ip = 0; ip < fsihist_.nvcvert; ip++) {
fsihist_.abspvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomLab;
fsihist_.abstpvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomNuc;
fsihist_.ipvert[ip] = nvect->FsiPartInfo(ip)->fPID;
fsihist_.iverti[ip] = nvect->FsiPartInfo(ip)->fVertStart;
fsihist_.ivertf[ip] = nvect->FsiPartInfo(ip)->fVertEnd;
fsihist_.dirvert[ip][0] = (float)nvect->FsiPartInfo(ip)->fDir.X();
fsihist_.dirvert[ip][1] = (float)nvect->FsiPartInfo(ip)->fDir.Y();
fsihist_.dirvert[ip][2] = (float)nvect->FsiPartInfo(ip)->fDir.Z();
}
fsihist_.fsiprob = nvect->Fsiprob;
}
neutcrscom_.crsx = nvect->Crsx;
neutcrscom_.crsy = nvect->Crsy;
neutcrscom_.crsz = nvect->Crsz;
neutcrscom_.crsphi = nvect->Crsphi;
neutcrscom_.crsq2 = nvect->Crsq2;
neuttarget_.numbndn = nvect->TargetA - nvect->TargetZ;
neuttarget_.numbndp = nvect->TargetZ;
neuttarget_.numfrep = nvect->TargetH;
neuttarget_.numatom = nvect->TargetA;
posinnuc_.ibound = nvect->Ibound;
// put empty nucleon FSI history (since it is not saved in the NeutVect format)
//Comment out as NEUT does not have the necessary proton FSI information yet
// nucleonfsihist_.nfnvert = 0;
// nucleonfsihist_.nfnstep = 0;
}
#endif
diff --git a/src/MINERvA/CMakeLists.txt b/src/MINERvA/CMakeLists.txt
index c8b8028..f74a9e9 100644
--- a/src/MINERvA/CMakeLists.txt
+++ b/src/MINERvA/CMakeLists.txt
@@ -1,131 +1,133 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
################################################################################
set(IMPLFILES
MINERvA_CCQE_XSec_1DQ2_antinu.cxx
MINERvA_CCQE_XSec_1DQ2_joint.cxx
MINERvA_CCQE_XSec_1DQ2_nu.cxx
MINERvA_CC0pi_XSec_1DEe_nue.cxx
MINERvA_CC0pi_XSec_1DQ2_nue.cxx
MINERvA_CC0pi_XSec_1DQ2_nu_proton.cxx
MINERvA_CC0pi_XSec_1DThetae_nue.cxx
MINERvA_CC1pi0_XSec_1DEnu_antinu.cxx
MINERvA_CC1pi0_XSec_1DQ2_antinu.cxx
MINERvA_CC1pi0_XSec_1Dpmu_antinu.cxx
MINERvA_CC1pi0_XSec_1Dppi0_antinu.cxx
MINERvA_CC1pi0_XSec_1DTpi0_antinu.cxx
MINERvA_CC1pi0_XSec_1Dth_antinu.cxx
MINERvA_CC1pi0_XSec_1Dthmu_antinu.cxx
MINERvA_CC1pip_XSec_1DTpi_20deg_nu.cxx
MINERvA_CC1pip_XSec_1DTpi_nu.cxx
MINERvA_CC1pip_XSec_1Dth_20deg_nu.cxx
MINERvA_CC1pip_XSec_1Dth_nu.cxx
MINERvA_CCNpip_XSec_1DEnu_nu.cxx
MINERvA_CCNpip_XSec_1DQ2_nu.cxx
MINERvA_CCNpip_XSec_1DTpi_nu.cxx
MINERvA_CCNpip_XSec_1Dpmu_nu.cxx
MINERvA_CCNpip_XSec_1Dth_nu.cxx
MINERvA_CCNpip_XSec_1Dthmu_nu.cxx
MINERvA_CCinc_XSec_2DEavq3_nu.cxx
MINERvA_CCinc_XSec_1Dx_ratio.cxx
MINERvA_CCinc_XSec_1DEnu_ratio.cxx
MINERvA_CCinc_XSec_1Dx_nu.cxx
MINERvA_CCinc_XSec_1DEnu_nu.cxx
MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.cxx
MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.cxx
+MINERvAUtils.cxx
MINERvA_SignalDef.cxx
)
set(HEADERFILES
MINERvA_CCQE_XSec_1DQ2_antinu.h
MINERvA_CCQE_XSec_1DQ2_joint.h
MINERvA_CCQE_XSec_1DQ2_nu.h
MINERvA_CC0pi_XSec_1DEe_nue.h
MINERvA_CC0pi_XSec_1DQ2_nue.h
MINERvA_CC0pi_XSec_1DQ2_nu_proton.h
MINERvA_CC0pi_XSec_1DThetae_nue.h
MINERvA_CC1pi0_XSec_1DEnu_antinu.h
MINERvA_CC1pi0_XSec_1DQ2_antinu.h
MINERvA_CC1pi0_XSec_1Dpmu_antinu.h
MINERvA_CC1pi0_XSec_1Dppi0_antinu.h
MINERvA_CC1pi0_XSec_1DTpi0_antinu.h
MINERvA_CC1pi0_XSec_1Dth_antinu.h
MINERvA_CC1pi0_XSec_1Dthmu_antinu.h
MINERvA_CC1pip_XSec_1DTpi_20deg_nu.h
MINERvA_CC1pip_XSec_1DTpi_nu.h
MINERvA_CC1pip_XSec_1Dth_20deg_nu.h
MINERvA_CC1pip_XSec_1Dth_nu.h
MINERvA_CCNpip_XSec_1DEnu_nu.h
MINERvA_CCNpip_XSec_1DQ2_nu.h
MINERvA_CCNpip_XSec_1DTpi_nu.h
MINERvA_CCNpip_XSec_1Dpmu_nu.h
MINERvA_CCNpip_XSec_1Dth_nu.h
MINERvA_CCNpip_XSec_1Dthmu_nu.h
MINERvA_CCinc_XSec_2DEavq3_nu.h
MINERvA_CCinc_XSec_1Dx_ratio.h
MINERvA_CCinc_XSec_1DEnu_ratio.h
MINERvA_CCinc_XSec_1Dx_nu.h
MINERvA_CCinc_XSec_1DEnu_nu.h
MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.h
MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.h
+MINERvAUtils.h
MINERvA_SignalDef.h
)
set(LIBNAME expMINERvA)
if(CMAKE_BUILD_TYPE MATCHES DEBUG)
add_library(${LIBNAME} STATIC ${IMPLFILES})
else(CMAKE_BUILD_TYPE MATCHES RELEASE)
add_library(${LIBNAME} SHARED ${IMPLFILES})
endif()
include_directories(${RWENGINE_INCLUDE_DIRECTORIES})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
include_directories(${CMAKE_SOURCE_DIR}/src/InputHandler)
include_directories(${CMAKE_SOURCE_DIR}/src/Genie)
include_directories(${CMAKE_SOURCE_DIR}/src/Utils)
include_directories(${CMAKE_SOURCE_DIR}/src/Reweight)
include_directories(${CMAKE_SOURCE_DIR}/src/Splines)
set_target_properties(${LIBNAME} PROPERTIES VERSION
"${ExtFit_VERSION_MAJOR}.${ExtFit_VERSION_MINOR}.${ExtFit_VERSION_REVISION}")
set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS})
if(DEFINED PROJECTWIDE_EXTRA_DEPENDENCIES)
add_dependencies(${LIBNAME} ${PROJECTWIDE_EXTRA_DEPENDENCIES})
endif()
install(TARGETS ${LIBNAME} DESTINATION lib)
#Can uncomment this to install the headers... but is it really neccessary?
#install(FILES ${HEADERFILES} DESTINATION include)
set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE)
diff --git a/src/MINERvA/MINERvAUtils.cxx b/src/MINERvA/MINERvAUtils.cxx
new file mode 100644
index 0000000..7670f45
--- /dev/null
+++ b/src/MINERvA/MINERvAUtils.cxx
@@ -0,0 +1,232 @@
+// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
+
+/*******************************************************************************
+* This file is part of NUISANCE.
+*
+* NUISANCE is free software: you can redistribute it and/or modify
+* it under the terms of the GNU General Public License as published by
+* the Free Software Foundation, either version 3 of the License, or
+* (at your option) any later version.
+*
+* NUISANCE is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
+*******************************************************************************/
+
+#include "MINERvAUtils.h"
+#include "FitParameters.h"
+#include "FitUtils.h"
+
+namespace FitPar{
+ double SciBarDensity = FitPar::Config().GetParD("SciBarDensity");
+ double SciBarRecoDist = FitPar::Config().GetParD("SciBarRecoDist");
+ double PenetratingMuonE = FitPar::Config().GetParD("PenetratingMuonEnergy");
+ double NumRangeSteps = FitPar::Config().GetParI("NumRangeSteps");
+}
+
+double MINERvAUtils::StoppedEfficiency(TH2D *effHist, FitParticle *nu, FitParticle *muon){
+
+ double eff = 0.;
+
+ if (!effHist) return eff;
+ eff = effHist->GetBinContent(effHist->FindBin(FitUtils::p(muon), FitUtils::th(nu, muon)/TMath::Pi()*180.));
+
+ return eff;
+}
+
+double MINERvAUtils::PenetratedEfficiency(FitParticle *nu, FitParticle *muon){
+
+ double eff = 0.;
+
+ if (FitUtils::th(nu, muon)/TMath::Pi()*180. > 50) eff = 0.;
+ if (FitUtils::p(muon) < 1.4) eff = 0.;
+
+ return eff;
+}
+
+double MINERvAUtils::BetheBlochCH(double E, double mass){
+
+ double beta2 = 1 - mass*mass/E/E;
+ double gamma = 1./sqrt(1-beta2);
+ double mass_ratio = PhysConst::mass_electron*1000./mass;
+ // Argh, have to remember to convert to MeV or you'll hate yourself!
+ double I2 = 68.7e-6*68.7e-6;
+
+ double w_max = 2*PhysConst::mass_electron*1000.*beta2*gamma*gamma;
+ w_max /= 1 + 2*gamma*mass_ratio + mass_ratio*mass_ratio;
+
+
+ // Values taken from the PDG for K = 0.307075 MeV mol-1 cm2, mean ionization energy I = 68.7 eV (Polystyrene)
+ // <Z/A> = 0.53768 (pdg.lbl.gov/AtomicNuclearProperties)
+ double log_term = log(2*PhysConst::mass_electron*1000.*beta2*gamma*gamma*w_max/I2);
+ double dedx = 0.307075*0.53768/beta2*(0.5*log_term - beta2);
+
+ return dedx;
+}
+
+
+// This function returns an estimate of the range of the particle in scintillator.
+// It uses crude integration and Bethe-Bloch to approximate the range.
+double MINERvAUtils::RangeInScintillator(FitParticle* particle, int nsteps){
+
+ // The particle energy
+ double E = particle->fP.E();
+ double M = particle->fP.M();
+ double Ek = E - M;
+
+ double step_size = Ek/float(nsteps+1);
+ double range = 0;
+
+ // Add an offset to make the integral a touch more accurate
+ Ek -= step_size/2.;
+ for (int i = 0; i < nsteps; ++i){
+
+ double dEdx = MINERvAUtils::BetheBlochCH(Ek+M, M);
+
+ Ek -= step_size;
+ // dEdx is -ve
+ range -= step_size/dEdx;
+ }
+
+ // Account for density of polystyrene
+ range /= FitPar::SciBarDensity;
+
+ // Range estimate is in cm
+ return range;
+}
+
+
+// Function to calculate the distance the particle travels in scintillator
+bool MINERvAUtils::PassesDistanceCut(FitParticle* beam, FitParticle* particle){
+
+ double dist = MINERvAUtils::RangeInScintillator(particle, FitPar::NumRangeSteps);
+ double zdist = dist*cos(FitUtils::th(beam, particle));
+
+ if (abs(zdist) < FitPar::SciBarRecoDist) return false;
+ return true;
+}
+
+
+// Function to return the MainTrk
+int MINERvAUtils::GetMainTrack(FitEvent *event, TH2D *effHist, FitParticle*& mainTrk, double& weight, bool penetrated){
+
+ FitParticle *nu = event->GetNeutrinoIn();
+ double highestMom = 0;
+ int index = 0;
+ mainTrk = NULL;
+
+ // Loop over particles
+ for (uint j = 2; j < event->Npart(); ++j){
+
+ // Final state only!
+ if (!(event->PartInfo(j))->fIsAlive) continue;
+ if (event->PartInfo(j)->fNEUTStatusCode != 0) continue;
+
+ int PID = event->PartInfo(j)->fPID;
+
+ // Only consider pions, muons for now
+ if (abs(PID) != 211 && abs(PID) != 13) continue;
+
+ // Ignore if higher momentum tracks available
+ if (FitUtils::p(event->PartInfo(j)) < highestMom) continue;
+
+ // Okay, now this is highest momentum
+ highestMom = FitUtils::p(event->PartInfo(j));
+ weight *= MINERvAUtils::StoppedEfficiency(effHist, nu, event->PartInfo(j));
+ index = j;
+ mainTrk = event->PartInfo(j);
+ } // end loop over particle stack
+
+ return index;
+}
+
+
+void MINERvAUtils::GetOtherTrackInfo(FitEvent *event, int mainIndex, int& nProtons, int& nPiMus, int& nVertex, FitParticle*& secondTrk){
+
+ // Reset everything
+ nPiMus = 0;
+ nProtons = 0;
+ nVertex = 0;
+ secondTrk = NULL;
+
+ double highestMom = 0.;
+
+ // Loop over particles
+ for (uint j = 2; j < event->Npart(); ++j){
+
+ // Don't re-count the main track
+ if (j == (uint)mainIndex) continue;
+
+ // Final state only!
+ if (!(event->PartInfo(j))->fIsAlive) continue;
+ if (event->PartInfo(j)->fNEUTStatusCode != 0) continue;
+
+ int PID = event->PartInfo(j)->fPID;
+
+ // Only consider pions, muons, protons
+ if (abs(PID) != 211 && PID != 2212 && abs(PID) != 13) continue;
+
+ // Must be reconstructed as a track in SciBooNE
+ if (MINERvAUtils::PassesDistanceCut(event->PartInfo(0), event->PartInfo(j))){
+
+ // Keep track of the second highest momentum track
+ if (FitUtils::p(event->PartInfo(j)) > highestMom){
+ highestMom = FitUtils::p(event->PartInfo(j));
+ secondTrk = event->PartInfo(j);
+ }
+
+ if (PID == 2212) nProtons += 1;
+ else nPiMus += 1;
+ } else nVertex += 1;
+
+ } // end loop over particle stack
+
+ return;
+}
+
+
+// NOTE: need to adapt this to allow for penetrating events...
+// Simpler, but gives the same results as in Hirade-san's thesis
+double MINERvAUtils::CalcThetaPr(FitEvent *event, FitParticle *main, FitParticle *second, bool penetrated){
+
+ FitParticle *nu = event->GetNeutrinoIn();
+
+ if (!main || !nu || !second) return -999;
+
+ // Construct the vector p_pr = (-p_mux, -p_muy, Enurec - pmucosthetamu)
+ // where p_mux, p_muy are the projections of the candidate muon momentum onto the x and y dimension respectively
+ double pmu = main->fP.Vect().Mag();
+ double pmu_x = main->fP.Vect().X();
+ double pmu_y = main->fP.Vect().Y();
+
+ if (penetrated){
+ pmu = 1400.;
+ double ratio = 1.4/main->fP.Vect().Mag();
+ TVector3 mod_mu = main->fP.Vect()*ratio;
+ pmu_x = mod_mu.X();
+ pmu_y = mod_mu.Y();
+ }
+
+ double Enuqe = FitUtils::EnuQErec(pmu/1000.,cos(FitUtils::th(nu, main)), 27., true)*1000.;
+ double p_pr_z = Enuqe - pmu*cos(FitUtils::th(nu, main));
+
+ TVector3 p_pr = TVector3(-pmu_x, -pmu_y, p_pr_z);
+ double thetapr = p_pr.Angle(second->fP.Vect())/TMath::Pi()*180.;
+
+ return thetapr;
+}
+
+double MINERvAUtils::CalcThetaPi(FitEvent *event, FitParticle *second){
+
+ FitParticle *nu = event->GetNeutrinoIn();
+
+ if (!second || !nu) return -999;
+
+ double thetapi = FitUtils::th(nu, second)/TMath::Pi()*180.;
+ return thetapi;
+}
+
diff --git a/src/MINERvA/MINERvAUtils.h b/src/MINERvA/MINERvAUtils.h
new file mode 100644
index 0000000..ed1976f
--- /dev/null
+++ b/src/MINERvA/MINERvAUtils.h
@@ -0,0 +1,71 @@
+// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
+
+/*******************************************************************************
+* This file is part of NUISANCE.
+*
+* NUISANCE is free software: you can redistribute it and/or modify
+* it under the terms of the GNU General Public License as published by
+* the Free Software Foundation, either version 3 of the License, or
+* (at your option) any later version.
+*
+* NUISANCE is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
+*******************************************************************************/
+
+#ifndef MINERvAUtils_H_SEEN
+#define MINERvAUtils_H_SEEN
+
+#include <math.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <cstring>
+#include <fstream>
+#include <iostream>
+#include <iostream>
+#include <numeric>
+#include <limits>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <TLorentzVector.h>
+#include <TH2D.h>
+
+#include "FitParticle.h"
+#include "FitEvent.h"
+#include "FitLogger.h"
+#include "StandardStacks.h"
+
+/*!
+ * \addtogroup Utils
+ * @{
+ */
+
+namespace FitPar {
+ extern double SciBarDensity;
+ extern double SciBarRecoDist;
+ extern double PenetratingMuonE;
+ extern double NumRangeSteps;
+}
+
+namespace MINERvAUtils {
+
+ double StoppedEfficiency(TH2D *effHist, FitParticle *nu, FitParticle *muon);
+ double PenetratedEfficiency(FitParticle *nu, FitParticle *muon);
+ double BetheBlochCH(double beta, double mass);
+ double RangeInScintillator(FitParticle* particle, int nsteps=50);
+
+ bool PassesDistanceCut(FitParticle* beam, FitParticle* particle);
+
+ int GetMainTrack(FitEvent *event, TH2D *effHist, FitParticle*& mainTrk, double& weight, bool penetrated=false);
+ void GetOtherTrackInfo(FitEvent *event, int mainIndex, int& nProtons, int& nPiMus, int& nVertex, FitParticle*& secondTrk);
+
+ double CalcThetaPr(FitEvent *event, FitParticle *main, FitParticle *second, bool penetrated=false);
+ double CalcThetaPi(FitEvent *event, FitParticle *second);
+
+}
+#endif
diff --git a/src/Reweight/GENIEWeightEngine.cxx b/src/Reweight/GENIEWeightEngine.cxx
index 3eef0e4..7483959 100644
--- a/src/Reweight/GENIEWeightEngine.cxx
+++ b/src/Reweight/GENIEWeightEngine.cxx
@@ -1,232 +1,237 @@
#include "GENIEWeightEngine.h"
GENIEWeightEngine::GENIEWeightEngine(std::string name) {
#ifdef __GENIE_ENABLED__
// Setup the NEUT Reweight engien
fCalcName = name;
LOG(FIT) << "Setting up GENIE RW : " << fCalcName << std::endl;
// Create RW Engine suppressing cout
StopTalking();
fGenieRW = new genie::rew::GReWeight();
// Get List of Vetos (Just for debugging)
std::string rw_engine_list = FitPar::Config().GetParS("FitWeight.fGenieRW_veto");
bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos;
bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos;
bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos;
bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos;
bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos;
bool xsec_resdec = rw_engine_list.find("hadro_res_decay") == std::string::npos;
bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos;
bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos;
bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos;
bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos;
bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos;
bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos;
bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos;
bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos;
bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos;
bool xsec_qeaxial = rw_engine_list.find("xsec_ccqe_axial") == std::string::npos;
// Now actually add the RW Calcs
if (xsec_ncel)
fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL);
if (xsec_ccqe)
fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE);
if (xsec_coh)
fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH);
if (xsec_nnres)
fGenieRW->AdoptWghtCalc("xsec_nonresbkg",
new genie::rew::GReWeightNonResonanceBkg);
if (xsec_nudis)
fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod);
if (xsec_resdec)
fGenieRW->AdoptWghtCalc("hadro_res_decay",
new genie::rew::GReWeightResonanceDecay);
if (xsec_fzone)
fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone);
if (xsec_intra)
fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke);
if (xsec_agky)
fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY);
if (xsec_qevec)
fGenieRW->AdoptWghtCalc("xsec_ccqe_vec",
new genie::rew::GReWeightNuXSecCCQEvec);
#if __GENIE_VERSION__ >= 212
if (xsec_qeaxial)
fGenieRW->AdoptWghtCalc("xsec_ccqe_axial",
new genie::rew::GReWeightNuXSecCCQEaxial);
#endif
if (xsec_dis)
fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS);
if (xsec_nc)
fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC);
if (xsec_ccres)
fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES);
if (xsec_ncres)
fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES);
if (xsec_nucqe)
fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM);
GReWeightNuXSecCCQE * rwccqe =
dynamic_cast<GReWeightNuXSecCCQE *> (fGenieRW->WghtCalc("xsec_ccqe"));
rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa);
// Default to include shape and normalization changes for CCRES (can be changed downstream if desired)
GReWeightNuXSecCCRES * rwccres =
dynamic_cast<GReWeightNuXSecCCRES *> (fGenieRW->WghtCalc("xsec_ccres"));
- rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv);
+ std::string marestype = FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode");
+ if (!marestype.compare("kModeNormAndMaMvShape")){ rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape); }
+ else if (!marestype.compare("kModeMaMv")){ rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv); }
+ else {
+ THROW("Unkown MARES Mode in GENIE Weight Engine : " << marestype );
+ }
// Default to include shape and normalization changes for NCRES (can be changed downstream if desired)
GReWeightNuXSecNCRES * rwncres =
dynamic_cast<GReWeightNuXSecNCRES *> (fGenieRW->WghtCalc("xsec_ncres"));
rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv);
// Default to include shape and normalization changes for DIS (can be changed downstream if desired)
GReWeightNuXSecDIS * rwdis =
dynamic_cast<GReWeightNuXSecDIS *> (fGenieRW->WghtCalc("xsec_dis"));
rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u);
// Final Reconfigure
fGenieRW->Reconfigure();
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
// allow cout again
StartTalking();
#else
ERR(FTL) << "GENIE RW NOT ENABLED" << std::endl;
#endif
};
void GENIEWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __GENIE_ENABLED__
// Get First enum
int nuisenum = Reweight::ConvDial(name, kGENIE);
// Setup Maps
fEnumIndex[nuisenum];// = std::vector<size_t>(0);
fNameIndex[name]; // = std::vector<size_t>(0);
// Split by commas
std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
for (uint i = 0; i < allnames.size(); i++) {
std::string singlename = allnames[i];
// Get RW
genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename);
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fGENIESysts.push_back(rwsyst);
// Initialize dial
std::cout << "Registering " << singlename << " from " << name << std::endl;
fGenieRW->Systematics().Init( fGENIESysts[index] );
// If Absolute
if (fIsAbsTwk) {
GSystUncertainty::Instance()->SetUncertainty( rwsyst, 1.0, 1.0 );
}
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
#endif
};
void GENIEWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __GENIE_ENABLED__
std::vector<size_t> indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
std::cout << "Setting GENIE Dial Value "
<< i << " " << indices[i] << " "
<< fGENIESysts[indices[i]] << std::endl;
fGenieRW->Systematics().Set( fGENIESysts[indices[i]], val);
}
#endif
}
void GENIEWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __GENIE_ENABLED__
std::vector<size_t> indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val);
}
#endif
}
void GENIEWeightEngine::Reconfigure(bool silent) {
#ifdef __GENIE_ENABLED__
// Hush now...
if (silent) StopTalking();
// Reconf
fGenieRW->Reconfigure();
fGenieRW->Print();
// Shout again
if (silent) StartTalking();
#endif
}
double GENIEWeightEngine::CalcWeight(BaseFitEvt* evt) {
double rw_weight = 1.0;
#ifdef __GENIE_ENABLED__
// Skip Non GENIE
if (evt->fType != kGENIE) return 1.0;
// Make nom weight
if (!evt) {
THROW("evt not found : " << evt);
}
if (!(evt->genie_event)) {
THROW("evt->genie_event not found!" << evt->genie_event);
}
if (!(evt->genie_event->event)) {
THROW("evt->genie_event->event GHepRecord not found!" << (evt->genie_event->event));
}
if (!fGenieRW) {
THROW("GENIE RW Not Found!" << fGenieRW);
}
rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event));
// std::cout << "Returning GENIE Weight for electron scattering = " << rw_weight << std::endl;
#endif
// Return rw_weight
return rw_weight;
}
diff --git a/src/Reweight/NUISANCESyst.cxx b/src/Reweight/NUISANCESyst.cxx
index c279f67..44366e1 100644
--- a/src/Reweight/NUISANCESyst.cxx
+++ b/src/Reweight/NUISANCESyst.cxx
@@ -1,51 +1,52 @@
#include "NUISANCESyst.h"
int Reweight::ConvertNUISANCEDial(std::string type) {
for (int i = kUnkownNUISANCEDial + 1; i < kNUISANCEDial_LAST; i++) {
if (!type.compare(ConvNUISANCEDial(i).c_str())) {
return i;
}
}
return kUnkownNUISANCEDial;
};
std::string Reweight::ConvNUISANCEDial(int type) {
switch (type) {
case kGaussianCorr_CCQE_norm: { return "GaussianCorr_CCQE_norm"; }
case kGaussianCorr_CCQE_tilt: { return "GaussianCorr_CCQE_tilt"; }
case kGaussianCorr_CCQE_Pq0: { return "GaussianCorr_CCQE_Pq0"; }
case kGaussianCorr_CCQE_Wq0: { return "GaussianCorr_CCQE_Wq0"; }
case kGaussianCorr_CCQE_Pq3: { return "GaussianCorr_CCQE_Pq3"; }
case kGaussianCorr_CCQE_Wq3: { return "GaussianCorr_CCQE_Wq3"; }
case kGaussianCorr_2p2h_norm: { return "GaussianCorr_2p2h_norm"; }
case kGaussianCorr_2p2h_tilt: { return "GaussianCorr_2p2h_tilt"; }
case kGaussianCorr_2p2h_Pq0: { return "GaussianCorr_2p2h_Pq0"; }
case kGaussianCorr_2p2h_Wq0: { return "GaussianCorr_2p2h_Wq0"; }
case kGaussianCorr_2p2h_Pq3: { return "GaussianCorr_2p2h_Pq3"; }
case kGaussianCorr_2p2h_Wq3: { return "GaussianCorr_2p2h_Wq3"; }
case kGaussianCorr_2p2h_PPandNN_norm: { return "GaussianCorr_2p2h_PPandNN_norm"; }
case kGaussianCorr_2p2h_PPandNN_tilt: { return "GaussianCorr_2p2h_PPandNN_tilt"; }
case kGaussianCorr_2p2h_PPandNN_Pq0: { return "GaussianCorr_2p2h_PPandNN_Pq0"; }
case kGaussianCorr_2p2h_PPandNN_Wq0: { return "GaussianCorr_2p2h_PPandNN_Wq0"; }
case kGaussianCorr_2p2h_PPandNN_Pq3: { return "GaussianCorr_2p2h_PPandNN_Pq3"; }
case kGaussianCorr_2p2h_PPandNN_Wq3: { return "GaussianCorr_2p2h_PPandNN_Wq3"; }
case kGaussianCorr_2p2h_NP_norm: { return "GaussianCorr_2p2h_NP_norm"; }
case kGaussianCorr_2p2h_NP_tilt: { return "GaussianCorr_2p2h_NP_tilt"; }
case kGaussianCorr_2p2h_NP_Pq0: { return "GaussianCorr_2p2h_NP_Pq0"; }
case kGaussianCorr_2p2h_NP_Wq0: { return "GaussianCorr_2p2h_NP_Wq0"; }
case kGaussianCorr_2p2h_NP_Pq3: { return "GaussianCorr_2p2h_NP_Pq3"; }
case kGaussianCorr_2p2h_NP_Wq3: { return "GaussianCorr_2p2h_NP_Wq3"; }
case kGaussianCorr_CC1pi_norm: { return "GaussianCorr_CC1pi_norm"; }
case kGaussianCorr_CC1pi_tilt: { return "GaussianCorr_CC1pi_tilt"; }
case kGaussianCorr_CC1pi_Pq0: { return "GaussianCorr_CC1pi_Pq0"; }
case kGaussianCorr_CC1pi_Wq0: { return "GaussianCorr_CC1pi_Wq0"; }
case kGaussianCorr_CC1pi_Pq3: { return "GaussianCorr_CC1pi_Pq3"; }
case kGaussianCorr_CC1pi_Wq3: { return "GaussianCorr_CC1pi_Wq3"; }
case kGaussianCorr_AllowSuppression: { return "GaussianCorr_AllowSuppression"; }
+ case kModeNorm_NormRES: { return "NormRES"; }
default: return "unknown_nuisance_dial";
}
};
diff --git a/src/Reweight/NUISANCESyst.h b/src/Reweight/NUISANCESyst.h
index 6315e23..296c956 100644
--- a/src/Reweight/NUISANCESyst.h
+++ b/src/Reweight/NUISANCESyst.h
@@ -1,48 +1,50 @@
#ifndef NUISANCESyst_H
#define NUISANCESyst_H
#include "GeneralUtils.h"
namespace Reweight {
enum NUISANCESyst {
kUnkownNUISANCEDial = 0,
kGaussianCorr_CCQE_norm,
kGaussianCorr_CCQE_tilt,
kGaussianCorr_CCQE_Pq0,
kGaussianCorr_CCQE_Wq0,
kGaussianCorr_CCQE_Pq3,
kGaussianCorr_CCQE_Wq3,
kGaussianCorr_2p2h_norm,
kGaussianCorr_2p2h_tilt,
kGaussianCorr_2p2h_Pq0,
kGaussianCorr_2p2h_Wq0,
kGaussianCorr_2p2h_Pq3,
kGaussianCorr_2p2h_Wq3,
kGaussianCorr_2p2h_PPandNN_norm,
kGaussianCorr_2p2h_PPandNN_tilt,
kGaussianCorr_2p2h_PPandNN_Pq0,
kGaussianCorr_2p2h_PPandNN_Wq0,
kGaussianCorr_2p2h_PPandNN_Pq3,
kGaussianCorr_2p2h_PPandNN_Wq3,
kGaussianCorr_2p2h_NP_norm,
kGaussianCorr_2p2h_NP_tilt,
kGaussianCorr_2p2h_NP_Pq0,
kGaussianCorr_2p2h_NP_Wq0,
kGaussianCorr_2p2h_NP_Pq3,
kGaussianCorr_2p2h_NP_Wq3,
kGaussianCorr_CC1pi_norm,
kGaussianCorr_CC1pi_tilt,
kGaussianCorr_CC1pi_Pq0,
kGaussianCorr_CC1pi_Wq0,
kGaussianCorr_CC1pi_Pq3,
kGaussianCorr_CC1pi_Wq3,
kGaussianCorr_AllowSuppression,
+ kModeNorm_NormRES,
+
kNUISANCEDial_LAST
};
int ConvertNUISANCEDial(std::string type);
std::string ConvNUISANCEDial(int type);
};
#endif
diff --git a/src/Reweight/NUISANCEWeightCalcs.cxx b/src/Reweight/NUISANCEWeightCalcs.cxx
index d3f36e9..1f3ddf2 100644
--- a/src/Reweight/NUISANCEWeightCalcs.cxx
+++ b/src/Reweight/NUISANCEWeightCalcs.cxx
@@ -1,276 +1,329 @@
#include "NUISANCEWeightCalcs.h"
+ModeNormCalc::ModeNormCalc(){
+ fNormRES = 1.0;
+}
+
+double ModeNormCalc::CalcWeight(BaseFitEvt* evt) {
+ FitEvent* fevt = static_cast<FitEvent*>(evt);
+ int mode = abs(fevt->Mode);
+ double w = 1.0;
+
+ if (mode == 11 or mode == 12 or mode == 13){
+ w *= fNormRES;
+ }
+
+ return w;
+}
+
+void ModeNormCalc::SetDialValue(std::string name, double val) {
+ SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
+}
+
+void ModeNormCalc::SetDialValue(int rwenum, double val) {
+
+ int curenum = rwenum % 1000;
+
+ // Check Handled
+ if (!IsHandled(curenum)) return;
+ if (curenum == kModeNorm_NormRES) fNormRES = val;
+}
+
+bool ModeNormCalc::IsHandled(int rwenum) {
+
+ int curenum = rwenum % 1000;
+ switch (curenum) {
+ case kModeNorm_NormRES: return true;
+ default: return false;
+ }
+}
+
+
+
+
+
+
+
GaussianModeCorr::GaussianModeCorr() {
// Init
fApply_CCQE = false;
fGausVal_CCQE[kPosNorm] = 0.0;
fGausVal_CCQE[kPosTilt] = 0.0;
fGausVal_CCQE[kPosPq0] = 1.0;
fGausVal_CCQE[kPosWq0] = 1.0;
fGausVal_CCQE[kPosPq3] = 1.0;
fGausVal_CCQE[kPosWq3] = 1.0;
fApply_2p2h = false;
fGausVal_2p2h[kPosNorm] = 0.0;
fGausVal_2p2h[kPosTilt] = 0.0;
fGausVal_2p2h[kPosPq0] = 1.0;
fGausVal_2p2h[kPosWq0] = 1.0;
fGausVal_2p2h[kPosPq3] = 1.0;
fGausVal_2p2h[kPosWq3] = 1.0;
fApply_2p2h_PPandNN = false;
fGausVal_2p2h_PPandNN[kPosNorm] = 0.0;
fGausVal_2p2h_PPandNN[kPosTilt] = 0.0;
fGausVal_2p2h_PPandNN[kPosPq0] = 1.0;
fGausVal_2p2h_PPandNN[kPosWq0] = 1.0;
fGausVal_2p2h_PPandNN[kPosPq3] = 1.0;
fGausVal_2p2h_PPandNN[kPosWq3] = 1.0;
fApply_2p2h_NP = false;
fGausVal_2p2h_NP[kPosNorm] = 0.0;
fGausVal_2p2h_NP[kPosTilt] = 0.0;
fGausVal_2p2h_NP[kPosPq0] = 1.0;
fGausVal_2p2h_NP[kPosWq0] = 1.0;
fGausVal_2p2h_NP[kPosPq3] = 1.0;
fGausVal_2p2h_NP[kPosWq3] = 1.0;
fApply_CC1pi = false;
fGausVal_CC1pi[kPosNorm] = 0.0;
fGausVal_CC1pi[kPosTilt] = 0.0;
fGausVal_CC1pi[kPosPq0] = 1.0;
fGausVal_CC1pi[kPosWq0] = 1.0;
fGausVal_CC1pi[kPosPq3] = 1.0;
fGausVal_CC1pi[kPosWq3] = 1.0;
fDebugStatements = FitPar::Config().GetParB("GaussianModeCorr_DEBUG");
}
double GaussianModeCorr::CalcWeight(BaseFitEvt* evt) {
FitEvent* fevt = static_cast<FitEvent*>(evt);
double rw_weight = 1.0;
// Get Neutrino
+ if (!fevt->Npart()){
+ THROW("NO particles found in stack!");
+ }
FitParticle* pnu = fevt->PartInfo(0);
+ if (!pnu){
+ THROW("NO Starting particle found in stack!");
+ }
int pdgnu = pnu->fPID;
FitParticle* plep = fevt->GetHMFSParticle(abs(pdgnu) - 1);
if (!plep) return 1.0;
TLorentzVector q = pnu->fP - plep->fP;
// Extra q0,q3
double q0 = fabs(q.E()) / 1.E3;
double q3 = fabs(q.Vect().Mag()) / 1.E3;
int initialstate = -1; // Undef
if (abs(fevt->Mode) == 2) {
int npr = 0;
int nne = 0;
for (UInt_t j = 0; j < fevt->Npart(); j++) {
if ((fevt->PartInfo(j))->fIsAlive) continue;
if (fevt->PartInfo(j)->fPID == 2212) npr++;
else if (fevt->PartInfo(j)->fPID == 2112) nne++;
}
// std::cout << "PN State = " << npr << " " << nne << std::endl;
if (fevt->Mode == 2 and npr == 1 and nne == 1) {
initialstate = 2;
} else if (fevt->Mode == 2 and ((npr == 0 and nne == 2) or (npr == 2 and nne == 0))) {
initialstate = 1;
}
}
// std::cout << "Got q0 q3 = " << q0 << " " << q3 << std::endl;
// Apply weighting
if (fApply_CCQE and abs(fevt->Mode) == 1) {
if (fDebugStatements) std::cout << "Getting CCQE Weight" << std::endl;
- rw_weight *= GetGausWeight(q0, q3, fGausVal_CCQE);
+ double g = GetGausWeight(q0, q3, fGausVal_CCQE);
+ if (g < 1.0) g = 1.0;
+ rw_weight *= g;
}
if (fApply_2p2h and abs(fevt->Mode) == 2) {
if (fDebugStatements) std::cout << "Getting 2p2h Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h);
}
if (fApply_2p2h_PPandNN and abs(fevt->Mode) == 2 and initialstate == 1) {
if (fDebugStatements) std::cout << "Getting 2p2h PPandNN Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_PPandNN);
}
if (fApply_2p2h_NP and abs(fevt->Mode) == 2 and initialstate == 2) {
if (fDebugStatements) std::cout << "Getting 2p2h NP Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_NP);
}
if (fApply_CC1pi and abs(fevt->Mode) >= 11 and abs(fevt->Mode) <= 13) {
if (fDebugStatements) std::cout << "Getting CC1pi Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_CC1pi);
}
// if (fDebugStatements) std::cout << "Returning Weight " << rw_weight << std::endl;
return rw_weight;
}
double GaussianModeCorr::GetGausWeight(double q0, double q3, double vals[]) {
// // CCQE Without Suppression
// double Norm = 4.82788679036;
// double Tilt = 2.3501416116;
// double Pq0 = 0.363964889702;
// double Wq0 = 0.133976806938;
// double Pq3 = 0.431769740224;
// double Wq3 = 0.207666663434;
// // Also add for CCQE at the end
// return (w > 1.0) ? w : 1.0;
// // 2p2h with suppression
// double Norm = 15.967;
// double Tilt = -0.455655;
// double Pq0 = 0.214598;
// double Wq0 = 0.0291061;
// double Pq3 = 0.480194;
// double Wq3 = 0.134588;
+
double Norm = vals[kPosNorm];
double Tilt = vals[kPosTilt];
double Pq0 = vals[kPosPq0];
double Wq0 = vals[kPosWq0];
double Pq3 = vals[kPosPq3];
double Wq3 = vals[kPosWq3];
double a = cos(Tilt) * cos(Tilt) / (2 * Wq0 * Wq0);
a += sin(Tilt) * sin(Tilt) / (2 * Wq3 * Wq3);
double b = - sin(2 * Tilt) / (4 * Wq0 * Wq0);
b += sin(2 * Tilt) / (4 * Wq3 * Wq3);
double c = sin(Tilt) * sin(Tilt) / (2 * Wq0 * Wq0);
c += cos(Tilt) * cos(Tilt) / (2 * Wq3 * Wq3);
double w = Norm;
w *= exp(-a * (q0 - Pq0) * (q0 - Pq0));
w *= exp(+2.0 * b * (q0 - Pq0) * (q3 - Pq3));
w *= exp(-c * (q3 - Pq3) * (q3 - Pq3));
if (fDebugStatements) {
std::cout << "Applied Tilt " << Tilt << " " << cos(Tilt) << " " << sin(Tilt) << std::endl;
std::cout << "abc = " << a << " " << b << " " << c << std::endl;
std::cout << "Returning " << Norm << " " << Pq0 << " " << Wq0 << " " << Pq3 << " " << Wq3 << " " << w << std::endl;
}
return w;
}
void GaussianModeCorr::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void GaussianModeCorr::SetDialValue(int rwenum, double val) {
int curenum = rwenum % 1000;
// Check Handled
if (!IsHandled(curenum)) return;
// CCQE Setting
for (int i = kGaussianCorr_CCQE_norm; i <= kGaussianCorr_CCQE_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_CCQE_norm;
fGausVal_CCQE[index] = val;
fApply_CCQE = true;
}
}
// 2p2h Setting
for (int i = kGaussianCorr_2p2h_norm; i <= kGaussianCorr_2p2h_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_norm;
fGausVal_2p2h[index] = val;
fApply_2p2h = true;
}
}
// 2p2h_PPandNN Setting
for (int i = kGaussianCorr_2p2h_PPandNN_norm; i <= kGaussianCorr_2p2h_PPandNN_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_PPandNN_norm;
fGausVal_2p2h_PPandNN[index] = val;
fApply_2p2h_PPandNN = true;
}
}
// 2p2h_NP Setting
for (int i = kGaussianCorr_2p2h_NP_norm; i <= kGaussianCorr_2p2h_NP_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_NP_norm;
fGausVal_2p2h_NP[index] = val;
fApply_2p2h_NP = true;
}
}
// CC1pi Setting
for (int i = kGaussianCorr_CC1pi_norm; i <= kGaussianCorr_CC1pi_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_CC1pi_norm;
fGausVal_CC1pi[index] = val;
fApply_CC1pi = true;
}
}
}
bool GaussianModeCorr::IsHandled(int rwenum) {
int curenum = rwenum % 1000;
switch (curenum) {
case kGaussianCorr_CCQE_norm:
case kGaussianCorr_CCQE_tilt:
case kGaussianCorr_CCQE_Pq0:
case kGaussianCorr_CCQE_Wq0:
case kGaussianCorr_CCQE_Pq3:
case kGaussianCorr_CCQE_Wq3:
case kGaussianCorr_2p2h_norm:
case kGaussianCorr_2p2h_tilt:
case kGaussianCorr_2p2h_Pq0:
case kGaussianCorr_2p2h_Wq0:
case kGaussianCorr_2p2h_Pq3:
case kGaussianCorr_2p2h_Wq3:
case kGaussianCorr_2p2h_PPandNN_norm:
case kGaussianCorr_2p2h_PPandNN_tilt:
case kGaussianCorr_2p2h_PPandNN_Pq0:
case kGaussianCorr_2p2h_PPandNN_Wq0:
case kGaussianCorr_2p2h_PPandNN_Pq3:
case kGaussianCorr_2p2h_PPandNN_Wq3:
case kGaussianCorr_2p2h_NP_norm:
case kGaussianCorr_2p2h_NP_tilt:
case kGaussianCorr_2p2h_NP_Pq0:
case kGaussianCorr_2p2h_NP_Wq0:
case kGaussianCorr_2p2h_NP_Pq3:
case kGaussianCorr_2p2h_NP_Wq3:
case kGaussianCorr_CC1pi_norm:
case kGaussianCorr_CC1pi_tilt:
case kGaussianCorr_CC1pi_Pq0:
case kGaussianCorr_CC1pi_Wq0:
case kGaussianCorr_CC1pi_Pq3:
case kGaussianCorr_CC1pi_Wq3:
return true;
default:
return false;
}
}
diff --git a/src/Reweight/NUISANCEWeightCalcs.h b/src/Reweight/NUISANCEWeightCalcs.h
index 4d66e4a..95fa64b 100644
--- a/src/Reweight/NUISANCEWeightCalcs.h
+++ b/src/Reweight/NUISANCEWeightCalcs.h
@@ -1,85 +1,100 @@
#ifndef NUISANCE_WEIGHT_CALCS
#define NUISANCE_WEIGHT_CALCS
#include "FitEvent.h"
#include "GeneralUtils.h"
#include "BaseFitEvt.h"
#include "WeightUtils.h"
#include "NUISANCESyst.h"
#include "FitParameters.h"
using namespace Reweight;
+
class NUISANCEWeightCalc {
public:
NUISANCEWeightCalc() {};
virtual ~NUISANCEWeightCalc() {};
virtual double CalcWeight(BaseFitEvt* evt){return 1.0;};
virtual void SetDialValue(std::string name, double val){};
virtual void SetDialValue(int rwenum, double val){};
virtual bool IsHandled(int rwenum){return false;};
virtual void Print(){};
std::map<std::string, int> fDialNameIndex;
std::map<int, int> fDialEnumIndex;
std::vector<double> fDialValues;
std::string fName;
};
+class ModeNormCalc : public NUISANCEWeightCalc {
+ public:
+ ModeNormCalc();
+ ~ModeNormCalc(){};
+
+ double CalcWeight(BaseFitEvt* evt);
+ void SetDialValue(std::string name, double val);
+ void SetDialValue(int rwenum, double val);
+ bool IsHandled(int rwenum);
+
+ double fNormRES;
+ };
+
+
class GaussianModeCorr : public NUISANCEWeightCalc {
public:
GaussianModeCorr();
~GaussianModeCorr(){};
double CalcWeight(BaseFitEvt* evt);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
bool IsHandled(int rwenum);
double GetGausWeight(double q0, double q3, double vals[]);
// 5 pars describe the Gaussain
// 0 norm.
// 1 q0 pos
// 2 q0 width
// 3 q3 pos
// 4 q3 width
static const int kPosNorm = 0;
static const int kPosTilt = 1;
static const int kPosPq0 = 2;
static const int kPosWq0 = 3;
static const int kPosPq3 = 4;
static const int kPosWq3 = 5;
bool fApply_CCQE;
double fGausVal_CCQE[6];
bool fApply_2p2h;
double fGausVal_2p2h[6];
bool fApply_2p2h_PPandNN;
double fGausVal_2p2h_PPandNN[6];
bool fApply_2p2h_NP;
double fGausVal_2p2h_NP[6];
bool fApply_CC1pi;
double fGausVal_CC1pi[6];
bool fAllowSuppression;
bool fDebugStatements;
};
#endif
diff --git a/src/Reweight/NUISANCEWeightEngine.cxx b/src/Reweight/NUISANCEWeightEngine.cxx
index 44b301a..e7344ca 100644
--- a/src/Reweight/NUISANCEWeightEngine.cxx
+++ b/src/Reweight/NUISANCEWeightEngine.cxx
@@ -1,129 +1,129 @@
#include "NUISANCEWeightEngine.h"
#include "NUISANCEWeightCalcs.h"
NUISANCEWeightEngine::NUISANCEWeightEngine(std::string name) {
// Setup the NUISANCE Reweight engine
fCalcName = name;
LOG(FIT) << "Setting up NUISANCE Custom RW : " << fCalcName << std::endl;
// Load in all Weight Calculations
fWeightCalculators.push_back( new GaussianModeCorr() );
-
+ fWeightCalculators.push_back( new ModeNormCalc() );
// Set Abs Twk Config
fIsAbsTwk = true;
};
void NUISANCEWeightEngine::IncludeDial(std::string name, double startval) {
// Get First enum
int nuisenum = Reweight::ConvDial(name, kCUSTOM);
// Setup Maps
fEnumIndex[nuisenum];// = std::vector<size_t>(0);
fNameIndex[name]; // = std::vector<size_t>(0);
// Split by commas
std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
for (uint i = 0; i < allnames.size(); i++) {
std::string singlename = allnames[i];
// Get RW
int singleenum = Reweight::ConvDial(singlename, kCUSTOM);
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fNUISANCEEnums.push_back(singleenum);
// Initialize dial
std::cout << "Registering " << singlename << " from " << name << std::endl;
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
};
void NUISANCEWeightEngine::SetDialValue(int nuisenum, double val) {
std::vector<size_t> indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
}
}
void NUISANCEWeightEngine::SetDialValue(std::string name, double val) {
std::vector<size_t> indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
}
}
void NUISANCEWeightEngine::Reconfigure(bool silent) {
for (size_t i = 0; i < fNUISANCEEnums.size(); i++) {
for (std::vector<NUISANCEWeightCalc*>::iterator calciter = fWeightCalculators.begin();
calciter != fWeightCalculators.end(); calciter++) {
NUISANCEWeightCalc* nuiscalc = static_cast<NUISANCEWeightCalc*>(*calciter);
if (nuiscalc->IsHandled(fNUISANCEEnums[i])) {
nuiscalc->SetDialValue(fNUISANCEEnums[i], fValues[i]);
}
}
}
}
double NUISANCEWeightEngine::CalcWeight(BaseFitEvt * evt) {
double rw_weight = 1.0;
// Cast as usable class
for (std::vector<NUISANCEWeightCalc*>::iterator iter = fWeightCalculators.begin();
iter != fWeightCalculators.end(); iter++) {
NUISANCEWeightCalc* nuiscalc = static_cast<NUISANCEWeightCalc*>(*iter);
rw_weight *= nuiscalc->CalcWeight(evt);
}
// Return rw_weight
return rw_weight;
}
diff --git a/src/Routines/SplineRoutines.cxx b/src/Routines/SplineRoutines.cxx
index 4cbe86a..8a24297 100755
--- a/src/Routines/SplineRoutines.cxx
+++ b/src/Routines/SplineRoutines.cxx
@@ -1,2593 +1,2593 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "SplineRoutines.h"
void SplineRoutines::Init() {
fStrategy = "SaveEvents";
fRoutines.clear();
fCardFile = "";
fSampleFCN = NULL;
fRW = NULL;
fAllowedRoutines = ("SaveEvents,TestEvents,SaveSplineEvents");
};
SplineRoutines::~SplineRoutines() {
};
SplineRoutines::SplineRoutines(int argc, char* argv[]) {
// Initialise Defaults
Init();
nuisconfig configuration = Config::Get();
// Default containers
std::string cardfile = "";
std::string maxevents = "-1";
int errorcount = 0;
int verbocount = 0;
std::vector<std::string> xmlcmds;
std::vector<std::string> configargs;
// Make easier to handle arguments.
std::vector<std::string> args = GeneralUtils::LoadCharToVectStr(argc, argv);
ParserUtils::ParseArgument(args, "-c", fCardFile, true);
ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false);
ParserUtils::ParseArgument(args, "-n", maxevents, false, false);
ParserUtils::ParseArgument(args, "-f", fStrategy, false, false);
ParserUtils::ParseArgument(args, "-i", xmlcmds);
ParserUtils::ParseArgument(args, "-q", configargs);
ParserUtils::ParseCounter(args, "e", errorcount);
ParserUtils::ParseCounter(args, "v", verbocount);
ParserUtils::CheckBadArguments(args);
// Add extra defaults if none given
if (fCardFile.empty() and xmlcmds.empty()) {
ERR(FTL) << "No input supplied!" << std::endl;
throw;
}
if (fOutputFile.empty() and !fCardFile.empty()) {
fOutputFile = fCardFile + ".root";
ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl;
} else if (fOutputFile.empty()) {
ERR(FTL) << "No output file or cardfile supplied!" << std::endl;
throw;
}
// Configuration Setup =============================
// Check no comp key is available
nuiskey fCompKey;
if (Config::Get().GetNodes("nuiscomp").empty()) {
fCompKey = Config::Get().CreateNode("nuiscomp");
} else {
fCompKey = Config::Get().GetNodes("nuiscomp")[0];
}
if (!fCardFile.empty()) fCompKey.AddS("cardfile", fCardFile);
fCompKey.AddS("outputfile", fOutputFile);
if (!fStrategy.empty()) fCompKey.AddS("strategy", fStrategy);
// Load XML Cardfile
configuration.LoadConfig( fCompKey.GetS("cardfile"), "");
// Add CMD XML Structs
for (size_t i = 0; i < xmlcmds.size(); i++) {
configuration.AddXMLLine(xmlcmds[i]);
}
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
if (maxevents.compare("-1")) {
std::cout << "[ NUISANCE ] : Overriding " << "MAXEVENTS=" + maxevents << std::endl;
configuration.OverrideConfig("MAXEVENTS=" + maxevents);
}
// Finish configuration XML
configuration.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml");
// Add Error Verbo Lines
verbocount += Config::Get().GetParI("VERBOSITY");
errorcount += Config::Get().GetParI("ERROR");
std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl;
std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl;
// FitPar::log_verb = verbocount;
SETVERBOSITY(verbocount);
// ERR_VERB(errorcount);
// Starting Setup
// ---------------------------
SetupRWEngine();
return;
};
/*
Setup Functions
*/
//*************************************
void SplineRoutines::SetupRWEngine() {
//*************************************
fRW = new FitWeight("splineweight");
// std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
// Add Parameters
for (size_t i = 0; i < parameterkeys.size(); i++) {
nuiskey key = parameterkeys[i];
std::string parname = key.GetS("name");
std::string partype = key.GetS("type");
double nom = key.GetD("nominal");
fRW->IncludeDial( key.GetS("name"),
FitBase::ConvDialType(key.GetS("type")), nom);
fRW->SetDialValue( key.GetS("name"), key.GetD("nominal") );
}
fRW->Reconfigure();
return;
}
/*
Fitting Functions
*/
//*************************************
void SplineRoutines::UpdateRWEngine(std::map<std::string, double>& updateVals) {
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string name = fParams[i];
if (updateVals.find(name) == updateVals.end()) continue;
fRW->SetDialValue(name, updateVals.at(name));
}
fRW->Reconfigure();
return;
}
//*************************************
void SplineRoutines::Run() {
//*************************************
std::cout << "Running " << std::endl;
// Parse given routines
fRoutines = GeneralUtils::ParseToStr(fStrategy, ",");
if (fRoutines.empty()) {
ERR(FTL) << "Trying to run ComparisonRoutines with no routines given!" << std::endl;
throw;
}
for (size_t i = 0; i < fRoutines.size(); i++) {
LOG(FIT) << "Running Routine: " << fRoutines[i] << std::endl;
std::string rout = fRoutines[i];
if (!rout.compare("SaveEvents")) SaveEvents();
else if (!rout.compare("TestEvents")) TestEvents();
else if (!rout.compare("GenerateEventSplines")) { GenerateEventWeights(); BuildEventSplines(); }
else if (!rout.compare("GenerateEventWeights")) { GenerateEventWeights(); }
else if (!rout.compare("GenerateEventWeightChunks")){ GenerateEventWeightChunks(FitPar::Config().GetParI("spline_procchunk")); }
else if (!rout.compare("BuildEventSplines")) { BuildEventSplines(); }
else if (!rout.compare("TestSplines_1DEventScan")) TestSplines_1DEventScan();
else if (!rout.compare("TestSplines_NDEventThrow")) TestSplines_NDEventThrow();
else if (!rout.compare("SaveSplinePlots")) SaveSplinePlots();
else if (!rout.compare("TestSplines_1DLikelihoodScan")) TestSplines_1DLikelihoodScan();
else if (!rout.compare("TestSplines_NDLikelihoodThrow")) TestSplines_NDLikelihoodThrow();
else if (!rout.compare("BuildEventSplinesChunks")) {
int chunk = FitPar::Config().GetParI("spline_procchunk");
BuildEventSplines(chunk);
} else if (!rout.compare("MergeEventSplinesChunks")) {
MergeEventSplinesChunks();
}
}
}
//*************************************
void SplineRoutines::SaveEvents() {
//*************************************
if (fRW) delete fRW;
SetupRWEngine();
fRW->Reconfigure();
fRW->Print();
// Generate a set of nominal events
// Method, Loop over inputs, create input handler, then create a ttree
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
for (size_t i = 0; i < eventkeys.size(); i++) {
nuiskey key = eventkeys.at(i);
// Get I/O
std::string inputfilename = key.GetS("input");
if (inputfilename.empty()) {
ERR(FTL) << "No input given for set of input events!" << std::endl;
throw;
}
std::string outputfilename = key.GetS("output");
if (outputfilename.empty()) {
outputfilename = inputfilename + ".nuisance.root";
ERR(FTL) << "No output give for set of output events! Saving to "
<< outputfilename << std::endl;
}
// Make new outputfile
TFile* outputfile = new TFile(outputfilename.c_str(), "RECREATE");
outputfile->cd();
// Make a new input handler
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(inputfilename, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inptype =
InputUtils::ParseInputType(file_descriptor[0]);
InputHandlerBase* input = InputUtils::CreateInputHandler("eventsaver", inptype, file_descriptor[1]);
// Get info from inputhandler
int nevents = input->GetNEvents();
int countwidth = (nevents / 10);
FitEvent* nuisevent = input->FirstNuisanceEvent();
// Setup a TTree to save the event
outputfile->cd();
TTree* eventtree = new TTree("nuisance_events", "nuisance_events");
nuisevent->AddBranchesToTree(eventtree);
// Loop over all events and fill the TTree
int icount = 0;
// int countwidth = nevents / 5;
while (nuisevent) {
// Get Event Weight
nuisevent->RWWeight = fRW->CalcWeight(nuisevent);
// if (nuisevent->RWWeight != 1.0){
// std::cout << "Weight = " << nuisevent->RWWeight << std::endl;
// }
// Save everything
eventtree->Fill();
// Logging
if (icount % countwidth == 0) {
LOG(REC) << "Saved " << icount << "/" << nevents
<< " nuisance events. [M, W] = ["
<< nuisevent->Mode << ", " << nuisevent->RWWeight << "]" << std::endl;
}
// iterate
nuisevent = input->NextNuisanceEvent();
- i++;
+ icount++;
}
// Save flux and close file
outputfile->cd();
eventtree->Write();
input->GetFluxHistogram()->Write("nuisance_fluxhist");
input->GetEventHistogram()->Write("nuisance_eventhist");
// Close Output
outputfile->Close();
// Delete Inputs
delete input;
}
// remove Keys
eventkeys.clear();
// Finished
LOG(FIT) << "Finished processing all nuisance events." << std::endl;
}
//*************************************
void SplineRoutines::TestEvents() {
//*************************************
LOG(FIT) << "Testing events." << std::endl;
// Create a new file for the test samples
if (!fOutputRootFile) {
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
}
// Loop over all tests
int count = 0;
std::vector<nuiskey> testkeys = Config::QueryKeys("sampletest");
for (std::vector<nuiskey>::iterator iter = testkeys.begin();
iter != testkeys.end(); iter++) {
nuiskey key = (*iter);
// 0. Create new measurement list
std::list<MeasurementBase*> samplelist;
// 1. Build Sample From Events
std::string samplename = key.GetS("name");
std::string eventsid = key.GetS("inputid");
nuiskey eventskey = Config::QueryLastKey("events", "id=" + eventsid);
std::string rawfile = eventskey.GetS("input");
LOG(FIT) << "Creating sample " << samplename << std::endl;
MeasurementBase* rawsample = SampleUtils::CreateSample(samplename, rawfile, "", "", FitBase::GetRW());
// 2. Build Sample From Nuisance Events
std::string eventsfile = eventskey.GetS("output");
LOG(FIT) << "Creating Fit Eevnt Sample " << samplename << " " << eventsfile << std::endl;
MeasurementBase* nuissample = SampleUtils::CreateSample(samplename, "FEVENT:" + eventsfile, "", "", FitBase::GetRW());
// 3. Make some folders to save stuff
TDirectory* sampledir = (TDirectory*) fOutputRootFile->mkdir(Form((samplename + "_test_%d").c_str(), count));
TDirectory* rawdir = (TDirectory*) sampledir->mkdir("raw");
TDirectory* nuisancedir = (TDirectory*) sampledir->mkdir("nuisance");
TDirectory* difdir = (TDirectory*) sampledir->mkdir("difference");
// 4. Reconfigure both
rawdir->cd();
rawsample->Reconfigure();
rawsample->Write();
nuisancedir->cd();
nuissample->Reconfigure();
nuissample->Write();
// 4. Compare Raw to Nuisance Events
// Loop over all keyse
TIter next(rawdir->GetListOfKeys());
TKey *dirkey;
while ((dirkey = (TKey*)next())) {
// If not a 1D/2D histogram skip
TClass *cl = gROOT->GetClass(dirkey->GetClassName());
if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue;
// Get TH1* from both dir
TH1 *rawplot = (TH1*)rawdir->Get(dirkey->GetName());
TH1 *nuisanceplot = (TH1*)nuisancedir->Get(dirkey->GetName());
// Take Difference
nuisanceplot->Add(rawplot, -1.0);
// Save to dif folder
difdir->cd();
nuisanceplot->Write();
}
// 5. Tidy Up
samplelist.clear();
// Iterator
count++;
}
}
void SplineRoutines::GenerateEventWeightChunks(int procchunk) {
if (fRW) delete fRW;
SetupRWEngine();
// Setup the spline reader
SplineWriter* splwrite = new SplineWriter(fRW);
std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
// Add splines to splinewriter
for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
iter != splinekeys.end(); iter++) {
nuiskey splkey = (*iter);
// Add Spline Info To Reader
splwrite->AddSpline(splkey);
}
splwrite->SetupSplineSet();
// Event Loop
// Loop over all events and calculate weights for each parameter set.
// Generate a set of nominal events
// Method, Loop over inputs, create input handler, then create a ttree
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
for (size_t i = 0; i < eventkeys.size(); i++) {
nuiskey key = eventkeys.at(i);
// Get I/O
std::string inputfilename = key.GetS("input");
if (inputfilename.empty()) {
ERR(FTL) << "No input given for set of input events!" << std::endl;
throw;
}
std::string outputfilename = key.GetS("output");
if (outputfilename.empty()) {
outputfilename = inputfilename + ".nuisance.root";
ERR(FTL) << "No output give for set of output events! Saving to "
<< outputfilename << std::endl;
}
outputfilename += ".weights.root";
// Make new outputfile
TFile* outputfile = new TFile(outputfilename.c_str(), "RECREATE");
outputfile->cd();
// Make a new input handler
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(inputfilename, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inptype =
InputUtils::ParseInputType(file_descriptor[0]);
InputHandlerBase* input = InputUtils::CreateInputHandler("eventsaver", inptype, file_descriptor[1]);
// Get info from inputhandler
int nevents = input->GetNEvents();
// int countwidth = (nevents / 1000);
FitEvent* nuisevent = input->FirstNuisanceEvent();
// Setup a TTree to save the event
outputfile->cd();
TTree* eventtree = new TTree("nuisance_events", "nuisance_events");
// Add a flag that allows just splines to be saved.
nuisevent->AddBranchesToTree(eventtree);
// Save the spline reader
splwrite->Write("spline_reader");
// Setup the spline TTree
TTree* weighttree = new TTree("weight_tree", "weight_tree");
splwrite->AddWeightsToTree(weighttree);
// Make container for all weights
int nweights = splwrite->GetNWeights();
// int npar = splwrite->GetNPars();
// double* weightcont = new double[nweights];
int lasttime = time(NULL);
// Load N Chunks of the Weights into Memory
// Split into N processing chunks
int nchunks = FitPar::Config().GetParI("spline_chunks");
if (nchunks <= 0) nchunks = 1;
if (nchunks >= nevents / 2) nchunks = nevents / 2;
std::cout << "Starting NChunks " << nchunks << std::endl;
for (int ichunk = 0; ichunk < nchunks; ichunk++) {
// Skip to only do one processing chunk
if (procchunk != -1 and procchunk != ichunk) continue;
LOG(FIT) << "On Processing Chunk " << ichunk << std::endl;
int neventsinchunk = nevents / nchunks;
int loweventinchunk = neventsinchunk * ichunk;
// int higheventinchunk = neventsinchunk * (ichunk + 1);
double** allweightcont = new double*[neventsinchunk];
for (int k = 0; k < neventsinchunk; k++){
allweightcont[k] = new double[nweights];
}
// Start Set Processing Here.
for (int iset = 0; iset < nweights; iset++) {
splwrite->ReconfigureSet(iset);
// Could reorder this to save the weightconts in order instead of reconfiguring per event.
// Loop over all events and fill the TTree
for (int i = 0; i < neventsinchunk; i++){
nuisevent = input->GetNuisanceEvent(i+loweventinchunk);
double w = splwrite->GetWeightForThisSet(nuisevent);
if (iset == 0){
allweightcont[i][0] = w;
} else {
allweightcont[i][iset] = w/allweightcont[i][0];
}
// Save everything
if (iset == 0) {
eventtree->Fill();
}
}
std::ostringstream timestring;
int timeelapsed = time(NULL) - lasttime;
if (timeelapsed) {
lasttime = time(NULL);
int setsleft = (nweights-iset-1) + (nweights * (nchunks - ichunk - 1));
float proj = (float(setsleft) *timeelapsed) / 60 / 60;
timestring << setsleft << " sets remaining. Last one took " << timeelapsed << ". " << proj << " hours remaining.";
}
LOG(REC) << "Processed Set " << iset << "/" << nweights <<" in chunk " << ichunk << "/" << nchunks << " " << timestring.str() << std::endl;
}
// Fill weights for this chunk into the TTree
for (int k = 0; k < neventsinchunk; k++){
splwrite->SetWeights(allweightcont[k]);
weighttree->Fill();
}
}
// at end of the chunk, when all sets have been done
// loop over the container and fill weights to ttree
outputfile->cd();
eventtree->Write();
weighttree->Write();
input->GetFluxHistogram()->Write("nuisance_fluxhist");
input->GetEventHistogram()->Write("nuisance_eventhist");
splwrite->Write("spline_reader");
outputfile->Close();
// Close Output
outputfile->Close();
// Delete Inputs
delete input;
}
// remove Keys
eventkeys.clear();
}
//*************************************
void SplineRoutines::GenerateEventWeights() {
//*************************************
if (fRW) delete fRW;
SetupRWEngine();
// Setup the spline reader
SplineWriter* splwrite = new SplineWriter(fRW);
std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
// Add splines to splinewriter
for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
iter != splinekeys.end(); iter++) {
nuiskey splkey = (*iter);
// Add Spline Info To Reader
splwrite->AddSpline(splkey);
}
splwrite->SetupSplineSet();
// Event Loop
// Loop over all events and calculate weights for each parameter set.
// Generate a set of nominal events
// Method, Loop over inputs, create input handler, then create a ttree
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
for (size_t i = 0; i < eventkeys.size(); i++) {
nuiskey key = eventkeys.at(i);
// Get I/O
std::string inputfilename = key.GetS("input");
if (inputfilename.empty()) {
ERR(FTL) << "No input given for set of input events!" << std::endl;
throw;
}
std::string outputfilename = key.GetS("output");
if (outputfilename.empty()) {
outputfilename = inputfilename + ".nuisance.root";
ERR(FTL) << "No output give for set of output events! Saving to "
<< outputfilename << std::endl;
}
outputfilename += ".weights.root";
// Make new outputfile
TFile* outputfile = new TFile(outputfilename.c_str(), "RECREATE");
outputfile->cd();
// Make a new input handler
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(inputfilename, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inptype =
InputUtils::ParseInputType(file_descriptor[0]);
InputHandlerBase* input = InputUtils::CreateInputHandler("eventsaver", inptype, file_descriptor[1]);
// Get info from inputhandler
int nevents = input->GetNEvents();
int countwidth = (nevents / 1000);
FitEvent* nuisevent = input->FirstNuisanceEvent();
// Setup a TTree to save the event
outputfile->cd();
TTree* eventtree = new TTree("nuisance_events", "nuisance_events");
// Add a flag that allows just splines to be saved.
nuisevent->AddBranchesToTree(eventtree);
// Save the spline reader
splwrite->Write("spline_reader");
// Setup the spline TTree
TTree* weighttree = new TTree("weight_tree", "weight_tree");
splwrite->AddWeightsToTree(weighttree);
// Make container for all weights
int nweights = splwrite->GetNWeights();
// int npar = splwrite->GetNPars();
double* weightcont = new double[nweights];
int lasttime = time(NULL);
// Could reorder this to save the weightconts in order instead of reconfiguring per event.
// Loop over all events and fill the TTree
while (nuisevent) {
// Calculate the weights for each parameter set
splwrite->GetWeightsForEvent(nuisevent, weightcont);
// Save everything
eventtree->Fill();
weighttree->Fill();
// Logging
if (i % countwidth == 0) {
std::ostringstream timestring;
int timeelapsed = time(NULL) - lasttime;
if (i != 0 and timeelapsed) {
lasttime = time(NULL);
int eventsleft = nevents - i;
float speed = float(countwidth) / float(timeelapsed);
float proj = (float(eventsleft) / float(speed)) / 60 / 60;
timestring << proj << " hours remaining.";
}
LOG(REC) << "Saved " << i << "/" << nevents << " nuisance spline weights. " << timestring.str() << std::endl;
}
// Iterate
i++;
nuisevent = input->NextNuisanceEvent();
}
// at end of the chunk, when all sets have been done
// loop over the container and fill weights to ttree
outputfile->cd();
eventtree->Write();
weighttree->Write();
input->GetFluxHistogram()->Write("nuisance_fluxhist");
input->GetEventHistogram()->Write("nuisance_eventhist");
splwrite->Write("spline_reader");
outputfile->Close();
// Close Output
outputfile->Close();
// Delete Inputs
delete input;
}
// remove Keys
eventkeys.clear();
}
//*************************************
void SplineRoutines::GenerateEventSplines() {
//*************************************
if (fRW) delete fRW;
SetupRWEngine();
// Setup the spline reader
SplineWriter* splwrite = new SplineWriter(fRW);
std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
// Add splines to splinewriter
for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
iter != splinekeys.end(); iter++) {
nuiskey splkey = (*iter);
// Add Spline Info To Reader
splwrite->AddSpline(splkey);
}
splwrite->SetupSplineSet();
// Make an ugly list for N cores
int ncores = FitPar::Config().GetParI("NCORES");//omp_get_max_threads();
std::vector<SplineWriter*> splwriterlist;
for (int i = 0; i < ncores; i++) {
SplineWriter* tmpwriter = new SplineWriter(fRW);
for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
iter != splinekeys.end(); iter++) {
nuiskey splkey = (*iter);
// Add Spline Info To Reader
tmpwriter->AddSpline(splkey);
}
tmpwriter->SetupSplineSet();
splwriterlist.push_back(tmpwriter);
}
// Event Loop
// Loop over all events and calculate weights for each parameter set.
// Generate a set of nominal events
// Method, Loop over inputs, create input handler, then create a ttree
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
for (size_t i = 0; i < eventkeys.size(); i++) {
nuiskey key = eventkeys.at(i);
// Get I/O
std::string inputfilename = key.GetS("input");
if (inputfilename.empty()) {
ERR(FTL) << "No input given for set of input events!" << std::endl;
throw;
}
std::string outputfilename = key.GetS("output");
if (outputfilename.empty()) {
outputfilename = inputfilename + ".nuisance.root";
ERR(FTL) << "No output give for set of output events! Saving to "
<< outputfilename << std::endl;
}
// Make new outputfile
TFile* outputfile = new TFile(outputfilename.c_str(), "RECREATE");
outputfile->cd();
// Make a new input handler
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(inputfilename, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inptype =
InputUtils::ParseInputType(file_descriptor[0]);
InputHandlerBase* input = InputUtils::CreateInputHandler("eventsaver", inptype, file_descriptor[1]);
// Get info from inputhandler
int nevents = input->GetNEvents();
int countwidth = (nevents / 1000);
FitEvent* nuisevent = input->FirstNuisanceEvent();
// Setup a TTree to save the event
outputfile->cd();
TTree* eventtree = new TTree("nuisance_events", "nuisance_events");
// Add a flag that allows just splines to be saved.
nuisevent->AddBranchesToTree(eventtree);
// Save the spline reader
splwrite->Write("spline_reader");
// Setup the spline TTree
TTree* weighttree = new TTree("weight_tree", "weight_tree");
splwrite->AddWeightsToTree(weighttree);
// Make container for all weights
int nweights = splwrite->GetNWeights();
double** weightcont = new double*[nevents];
for (int k = 0; k < nevents; k++) {
weightcont[k] = new double[nweights];
}
int npar = splwrite->GetNPars();
int lasttime = time(NULL);
// Could reorder this to save the weightconts in order instead of reconfiguring per event.
// Loop over all events and fill the TTree
while (nuisevent) {
// std::cout << "Fitting event " << i << std::endl;
// Calculate the weights for each parameter set
// splwrite->FitSplinesForEvent(nuisevent);
splwrite->GetWeightsForEvent(nuisevent, weightcont[i]);
bool hasresponse = false;
for (int j = 0; j < nweights; j++) {
if (weightcont[i][j] != 1.0) {
// std::cout << "Non Zero Weight at " << i << " " << j << std::endl;
hasresponse = true;
} else {
// std::cout << "Empty Weight at " << i << " " << j << std::endl;
}
}
if (!hasresponse) {
// std::cout << "Deleting flat response " << nuisevent->Mode << std::endl;
delete weightcont[i];
weightcont[i] = NULL;
}
// Save everything
eventtree->Fill();
weighttree->Fill();
// splinetree->Fill();
// nuisevent->Print();
// std::cout << "Done with event " << i << std::endl;
// Push weight sets into a the array
// sleep(4);
// Logging
if (i % countwidth == 0) {
std::ostringstream timestring;
int timeelapsed = time(NULL) - lasttime;
if (i != 0 and timeelapsed) {
lasttime = time(NULL);
int eventsleft = nevents - i;
float speed = float(countwidth) / float(timeelapsed);
float proj = (float(eventsleft) / float(speed)) / 60 / 60;
timestring << proj << " hours remaining.";
}
LOG(REC) << "Saved " << i << "/" << nevents << " nuisance spline weights. " << timestring.str() << std::endl;
}
// Iterate
i++;
nuisevent = input->NextNuisanceEvent();
}
outputfile->cd();
eventtree->Write();
weighttree->Write();
input->GetFluxHistogram()->Write("nuisance_fluxhist");
input->GetEventHistogram()->Write("nuisance_eventhist");
outputfile->Close();
outputfile = new TFile(outputfilename.c_str(), "UPDATE");
outputfile->cd();
weighttree = (TTree*) outputfile->Get("weight_tree");
// splwrite->ReadWeightsFromTree(weighttree);
// Divide weights container into Ncores.
// Parrallelise this loop checking for what core we are on.
// for (int i = 0; i < nevents; i++){
// splwriterlist[int(i / (nevents/4))]->FitSplinesForEvent(coeff);
// }
// // Now loop over weights tree
// for (int i = 0; i < weighttree->GetEntries(); i++) {
// weighttree->GetEntry(i);
// splwrite->FitSplinesForEvent();
// splinetree->Fill();
// if (i % countwidth == 0) {
// std::ostringstream timestring;
// int timeelapsed = time(NULL) - lasttime;
// if (i != 0 and timeelapsed) {
// lasttime = time(NULL);
// int eventsleft = nevents - i;
// float speed = float(countwidth) / float(timeelapsed);
// float proj = (float(eventsleft) / float(speed)) / 60 / 60;
// timestring << proj << " hours remaining.";
// }
// LOG(REC) << "Built " << i << "/" << nevents << " nuisance spline events. " << timestring.str() << std::endl;
// }
// }
// Get Splines
float** allcoeff = new float*[nevents];
for (int k = 0; k < nevents; k++) {
allcoeff[k] = new float[npar];
}
// #pragma omp parallel for num_threads(ncores)
for (int i = 0; i < nevents; i++) {
//#pragma omp atomic
// printf("Using Thread %d to build event %d \n", int(omp_get_thread_num()), (int)i );
// std::cout<< " -> Writer = " << splwriterlist[ i / (nevents/ncores) ] << std::endl;
// #pragma omp atomic
if (weightcont[i]) {
splwriterlist[ int(omp_get_thread_num()) ]->FitSplinesForEvent(weightcont[i], allcoeff[i]);
} else {
for (int j = 0; j < npar; j++) {
allcoeff[i][j] = float(0.0);
}
}
// splwrite->FitSplinesForEvent(weightcont[i], allcoeff[i]);
if (i % 500 == 0) {
if (LOG_LEVEL(REC)) {
printf("Using Thread %d to build event %d \n", int(omp_get_thread_num()), (int)i );
}
}
/*
std::ostringstream timestring;
int timeelapsed = time(NULL) - lasttime;
if (i != 0 and timeelapsed) {
lasttime = time(NULL);
int eventsleft = nevents - i;
float speed = float(countwidth) / float(timeelapsed);
float proj = (float(eventsleft) / float(speed)) / 60 / 60;
timestring << proj << " hours remaining.";
timestring << " Using Writer at " << i / (nevents/ncores) << " = " << splwriterlist[ i / (nevents/ncores) ] << std::endl;
}
LOG(REC) << "Built " << i << "/" << nevents << " nuisance spline events. " << timestring.str() << std::endl;
}
*/
}
// Save Splines into TTree
float* coeff = new float[npar];
outputfile->cd();
TTree* splinetree = new TTree("spline_tree", "spline_tree");
splinetree->Branch("SplineCoeff", coeff, Form("SplineCoeff[%d]/F", npar));
std::cout << "Saving to the allcoeff" << std::endl;
for (int k = 0; k < nevents; k++) {
for (int l = 0; l < npar; l++) {
coeff[l] = allcoeff[k][l];
}
std::cout << "Coeff 0, 1, 2 = " << coeff[0] << " " << coeff[1] << " " << coeff[2] << std::endl;
splinetree->Fill();
}
// Save flux and close file
outputfile->cd();
splinetree->Write();
// Delete the container.
for (int k = 0; k < nevents; k++) {
delete weightcont[k];
}
delete weightcont;
delete coeff;
// Close Output
outputfile->Close();
// Delete Inputs
delete input;
}
// remove Keys
eventkeys.clear();
}
//*************************************
void SplineRoutines::BuildEventSplines(int procchunk) {
//*************************************
if (fRW) delete fRW;
SetupRWEngine();
// Setup the spline reader
SplineWriter* splwrite = new SplineWriter(fRW);
std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
// Add splines to splinewriter
for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
iter != splinekeys.end(); iter++) {
nuiskey splkey = (*iter);
// Add Spline Info To Reader
splwrite->AddSpline(splkey);
}
splwrite->SetupSplineSet();
// Make an ugly list for N cores
int ncores = FitPar::Config().GetParI("spline_cores");//omp_get_max_threads();
if (ncores > omp_get_max_threads()) ncores = omp_get_max_threads();
if (ncores <= 0) ncores = 1;
std::vector<SplineWriter*> splwriterlist;
for (int i = 0; i < ncores; i++) {
SplineWriter* tmpwriter = new SplineWriter(fRW);
for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
iter != splinekeys.end(); iter++) {
nuiskey splkey = (*iter);
// Add Spline Info To Reader
tmpwriter->AddSpline(splkey);
}
tmpwriter->SetupSplineSet();
splwriterlist.push_back(tmpwriter);
}
// Event Loop
// Loop over all events and calculate weights for each parameter set.
// Generate a set of nominal events
// Method, Loop over inputs, create input handler, then create a ttree
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
for (size_t i = 0; i < eventkeys.size(); i++) {
nuiskey key = eventkeys.at(i);
// Get I/O
std::string inputfilename = key.GetS("input");
if (inputfilename.empty()) {
ERR(FTL) << "No input given for set of input events!" << std::endl;
throw;
}
std::string outputfilename = key.GetS("output");
if (outputfilename.empty()) {
outputfilename = inputfilename + ".nuisance.root";
ERR(FTL) << "No output give for set of output events! Saving to "
<< outputfilename << std::endl;
}
// Make new outputfile
TFile* outputfile;
if (procchunk == -1) outputfile = new TFile(outputfilename.c_str(), "RECREATE");
else outputfile = new TFile((outputfilename + std::string(Form(".coeffchunk_%d.root", procchunk))).c_str(), "RECREATE");
outputfile->cd();
// Get Weights File
TFile* weightsfile = new TFile((outputfilename + ".weights.root").c_str(), "READ");
TTree* weighttree = (TTree*) weightsfile->Get("weight_tree");
// Get SPLWRite Info
//splwrite->ReadWeightsFromTree(weighttree);
int nevents = weighttree->GetEntries();
// int countwidth = (nevents / 1000);
int nweights = splwrite->GetNWeights();
int npar = splwrite->GetNPars();
// Access Weights
double* eventweights = new double[nweights];
weighttree->SetBranchAddress("SplineWeights", eventweights);
// Make counter
// int lasttime = time(NULL);
// Setup Splines To Be Saved into TTree
outputfile->cd();
TTree* splinetree = new TTree("spline_tree", "spline_tree");
float* coeff = new float[npar];
splinetree->Branch("SplineCoeff", coeff, Form("SplineCoeff[%d]/F", npar));
// Load N Chunks of the Weights into Memory
// Split into N processing chunks
int nchunks = FitPar::Config().GetParI("spline_chunks");
if (nchunks <= 0) nchunks = 1;
if (nchunks >= nevents / 2) nchunks = nevents / 2;
std::cout << "Starting NChunks " << nchunks << std::endl;
sleep(1);
for (int ichunk = 0; ichunk < nchunks; ichunk++) {
// Skip to only do one processing chunk
if (procchunk != -1 and procchunk != ichunk) continue;
LOG(FIT) << "On Processing Chunk " << ichunk << std::endl;
int neventsinchunk = nevents / nchunks;
int loweventinchunk = neventsinchunk * ichunk;
// int higheventinchunk = neventsinchunk * (ichunk + 1);
// Build Chunk Containers for Event Weights
double** weightcont = new double*[nevents];
float** allcoeff = new float*[nevents];
// Load Chunks into Containers
for (int k = 0; k < neventsinchunk; k++) {
weighttree->GetEntry(loweventinchunk + k);
weightcont[k] = new double[nweights];
allcoeff[k] = new float[npar];
bool hasresponse = false;
for (int j = 0; j < nweights; j++) {
weightcont[k][j] = eventweights[j];
if (eventweights[j] != 1.0) hasresponse = true;
}
if (!hasresponse) delete weightcont[k];
}
// Loop over ncores and process chunks
// #pragma omp parallel for num_threads(ncores)
for (int k = 0; k < neventsinchunk; k++) {
if (weightcont[k]) {
splwriterlist[ int(omp_get_thread_num()) ]->FitSplinesForEvent(weightcont[k], allcoeff[k]);
} else {
for (int j = 0; j < npar; j++) {
allcoeff[k][j] = float(0.0);
}
}
if (k + loweventinchunk % 500 == 0) {
if (LOG_LEVEL(REC)) {
printf("Using Thread %d to build event %d in chunk %d \n", int(omp_get_thread_num()), (int) loweventinchunk + k, ichunk );
}
}
}
// Save Coeff To Tree
std::cout << "Saving coeffs to Tree in Chunk " << ichunk << std::endl;
for (int k = 0; k < neventsinchunk; k++) {
for (int l = 0; l < npar; l++) {
coeff[l] = allcoeff[k][l];
}
// std::cout << "Coeff 0, 1, 2 = " << coeff[0] << " " << coeff[1] << " " << coeff[2] << std::endl;
splinetree->Fill();
}
// Delete the container.
for (int k = 0; k < neventsinchunk; k++) {
if (weightcont[k]) delete weightcont[k];
if (allcoeff[k]) delete allcoeff[k];
}
delete allcoeff;
delete weightcont;
}
// Save flux and close file
outputfile->cd();
splinetree->Write();
if (procchunk == -1 or procchunk == 0) {
outputfile->cd();
splwrite->Write("spline_reader");
TTree* nuisanceevents = (TTree*) weightsfile->Get("nuisance_events");
nuisanceevents->CloneTree()->Write();
weighttree->CloneTree()->Write();
TH1D* nuisance_fluxhist = (TH1D*) weightsfile->Get("nuisance_fluxhist");
TH1D* nuisance_eventhist = (TH1D*) weightsfile->Get("nuisance_eventhist");
nuisance_fluxhist->Write("nuisance_fluxhist");
nuisance_eventhist->Write("nuisance_eventhist");
}
weightsfile->Close();
// Add option to build seperate chunks
// Close Output
outputfile->Close();
}
// remove Keys
eventkeys.clear();
}
void SplineRoutines::MergeEventSplinesChunks() {
if (fRW) delete fRW;
SetupRWEngine();
// Setup the spline reader
SplineWriter* splwrite = new SplineWriter(fRW);
std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
// Add splines to splinewriter
for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
iter != splinekeys.end(); iter++) {
nuiskey splkey = (*iter);
// Add Spline Info To Reader
splwrite->AddSpline(splkey);
}
splwrite->SetupSplineSet();
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
for (size_t i = 0; i < eventkeys.size(); i++) {
nuiskey key = eventkeys.at(i);
// Get I/O
std::string inputfilename = key.GetS("input");
if (inputfilename.empty()) {
ERR(FTL) << "No input given for set of input events!" << std::endl;
throw;
}
std::string outputfilename = key.GetS("output");
if (outputfilename.empty()) {
outputfilename = inputfilename + ".nuisance.root";
ERR(FTL) << "No output give for set of output events! Saving to "
<< outputfilename << std::endl;
}
// Make new outputfile
TFile* outputfile = new TFile(outputfilename.c_str(), "RECREATE");
outputfile->cd();
// Get Weights File
TFile* weightsfile = new TFile((outputfilename + ".weights.root").c_str(), "READ");
TTree* weighttree = (TTree*) weightsfile->Get("weight_tree");
// Get SPLWRite Info
//splwrite->ReadWeightsFromTree(weighttree);
int nevents = weighttree->GetEntries();
// int countwidth = (nevents / 1000);
// int nweights = splwrite->GetNWeights();
int npar = splwrite->GetNPars();
// Make counter
// int lasttime = time(NULL);
// Setup Splines To Be Saved into TTree
outputfile->cd();
TTree* splinetree = new TTree("spline_tree", "spline_tree");
float* coeff = new float[npar];
splinetree->Branch("SplineCoeff", coeff, Form("SplineCoeff[%d]/F", npar));
// Load N Chunks of the Weights into Memory
// Split into N processing chunks
int nchunks = FitPar::Config().GetParI("spline_chunks");
if (nchunks <= 0) nchunks = 1;
if (nchunks >= nevents / 2) nchunks = nevents / 2;
int neventsinchunk = nevents / nchunks;
for (int ichunk = 0; ichunk < nchunks; ichunk++) {
// Get Output File
TFile* chunkfile = new TFile( (outputfilename + std::string(Form(".coeffchunk_%d.root", ichunk))).c_str() );
// Get TTree for spline coeffchunk
TTree* splinetreechunk = (TTree*) chunkfile->Get("spline_tree");
// Set Branch Address to coeffchunk
float* coeffchunk = new float[npar];
splinetreechunk->SetBranchAddress("SplineCoeff", coeffchunk);
// Loop over nevents in chunk
for (int k = 0; k < neventsinchunk; k++) {
splinetreechunk->GetEntry(k);
for (int j = 0; j < npar; j++) {
coeff[j] = coeffchunk[j];
}
splinetree->Fill();
}
// Close up
chunkfile->Close();
delete coeffchunk;
std::cout << "Merged chunk " << ichunk << std::endl;
}
// Save flux and close file
outputfile->cd();
splinetree->Write();
outputfile->cd();
splwrite->Write("spline_reader");
TTree* nuisanceevents = (TTree*) weightsfile->Get("nuisance_events");
nuisanceevents->CloneTree()->Write();
weighttree->CloneTree()->Write();
TH1D* nuisance_fluxhist = (TH1D*) weightsfile->Get("nuisance_fluxhist");
TH1D* nuisance_eventhist = (TH1D*) weightsfile->Get("nuisance_eventhist");
nuisance_fluxhist->Write("nuisance_fluxhist");
nuisance_eventhist->Write("nuisance_eventhist");
weightsfile->Close();
// Add option to build seperate chunks
// Close Output
outputfile->Close();
}
// remove Keys
eventkeys.clear();
}
// void SplineRoutines::BuildSplineChunk(){
//}
// void SplineRoutines::MergeSplineChunks(){
//}
//*************************************
void SplineRoutines::MergeSplines() {
//*************************************
// Loop over all 'splinemerge' keys.
// Add them to the Merger.
// Call setup splines.
// Get the key with eventinput
// - remaining keys should have splineinput
// - Loop over number of entries.
// - FillEntry in merger.
// - Fill NUISANCEEvent into a new TTree.
SplineMerger* splmerge = new SplineMerger();
std::vector<nuiskey> splinekeys = Config::QueryKeys("splinemerge");
for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
iter != splinekeys.end(); iter++) {
nuiskey splkey = (*iter);
TFile* infile = new TFile(splkey.GetS("input").c_str(), "READ");
splmerge->AddSplineSetFromFile(infile);
}
splmerge->SetupSplineSet();
// Now get Event File
std::vector<nuiskey> eventkeys = Config::QueryKeys("eventmerge");
nuiskey key = eventkeys[0];
std::string inputfilename = key.GetS("input");
// Make a new input handler
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(inputfilename, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inptype =
InputUtils::ParseInputType(file_descriptor[0]);
InputHandlerBase* input = InputUtils::CreateInputHandler("eventsaver", inptype, file_descriptor[1]);
std::string outputfilename = key.GetS("output");
if (outputfilename.empty()) {
outputfilename = inputfilename + ".nuisance.root";
ERR(FTL) << "No output give for set of output events! Saving to "
<< outputfilename << std::endl;
}
// Make new outputfile
TFile* outputfile = new TFile(outputfilename.c_str(), "RECREATE");
outputfile->cd();
// Get info from inputhandler
int nevents = input->GetNEvents();
int countwidth = (nevents / 1000);
FitEvent* nuisevent = input->FirstNuisanceEvent();
// Setup a TTree to save the event
outputfile->cd();
TTree* eventtree = new TTree("nuisance_events", "nuisance_events");
// Add a flag that allows just splines to be saved.
nuisevent->AddBranchesToTree(eventtree);
// Save the spline reader
splmerge->Write("spline_reader");
// Setup the spline TTree
TTree* splinetree = new TTree("spline_tree", "spline_tree");
splmerge->AddCoefficientsToTree(splinetree);
int lasttime = time(NULL);
int i = 0;
// Loop over all events and fill the TTree
while (nuisevent) {
// Calculate the weights for each parameter set
splmerge->FillMergedSplines(i);
// Save everything
eventtree->Fill();
splinetree->Fill();
// Logging
if (i % countwidth == 0) {
std::ostringstream timestring;
int timeelapsed = time(NULL) - lasttime;
if (i != 0 and timeelapsed) {
lasttime = time(NULL);
int eventsleft = nevents - i;
float speed = float(countwidth) / float(timeelapsed);
float proj = (float(eventsleft) / float(speed)) / 60 / 60;
timestring << proj << " hours remaining.";
}
LOG(REC) << "Saved " << i << "/" << nevents << " nuisance spline events. " << timestring.str() << std::endl;
}
// Iterate
i++;
nuisevent = input->NextNuisanceEvent();
}
// Save flux and close file
outputfile->cd();
eventtree->Write();
splinetree->Write();
input->GetFluxHistogram()->Write("nuisance_fluxhist");
input->GetEventHistogram()->Write("nuisance_eventhist");
// Close Output
outputfile->Close();
// Delete Inputs
delete input;
}
//*************************************
void SplineRoutines::TestSplines_1DEventScan() {
//*************************************
// Setup RW Engine
if (fRW) delete fRW;
SetupRWEngine();
// Make a spline RW Engine too.
FitWeight* splweight = new FitWeight("splinerwaweight");
// std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
TH1D* parhisttemplate = new TH1D("parhist", "parhist", parameterkeys.size(), 0.0, float(parameterkeys.size()));
// Add Parameters
for (size_t i = 0; i < parameterkeys.size(); i++) {
nuiskey key = parameterkeys[i];
std::string parname = key.GetS("name");
std::string partype = key.GetS("type");
double nom = key.GetD("nominal");
parhisttemplate->SetBinContent(i + 1, nom);
parhisttemplate->GetXaxis()->SetBinLabel(i + 1, parname.c_str());
splweight->IncludeDial( key.GetS("name"),
kSPLINEPARAMETER, nom);
splweight->SetDialValue( key.GetS("name"), key.GetD("nominal") );
}
splweight->Reconfigure();
// Make a high resolution spline set.
std::vector<double> nomvals = fRW->GetDialValues();
// int testres = FitPar::Config().GetParI("spline_test_resolution");
std::vector< std::vector<double> > scanparset_vals;
std::vector< TH1D* > scanparset_hists;
// Loop over all params
// Add Parameters
for (size_t i = 0; i < parameterkeys.size(); i++) {
nuiskey key = parameterkeys[i];
// Get Par Name
std::string name = key.GetS("name");
if (!key.Has("low") or !key.Has("high") or !key.Has("step")) {
continue;
}
// Push Back Scan
double low = key.GetD("low");
double high = key.GetD("high");
double cur = low;
double step = key.GetD("step");
while (cur <= high) {
// Make new set
std::vector<double> newvals = nomvals;
newvals[i] = cur;
// Add to vects
scanparset_vals.push_back(newvals);
TH1D* parhist = (TH1D*)parhisttemplate->Clone();
for (size_t j = 0; j < newvals.size(); j++) {
parhist->SetBinContent(j + 1, newvals[j]);
}
scanparset_hists.push_back(parhist);
// Move to next one
cur += step;
}
}
// Print out the parameter set to test
for (uint i = 0; i < scanparset_vals.size(); i++) {
std::cout << "Parset " << i;
for (uint j = 0 ; j < scanparset_vals[i].size(); j++) {
std::cout << " " << scanparset_vals[i][j];
}
std::cout << std::endl;
}
// Weight holders
double* rawweights = new double[scanparset_vals.size()];
double* splweights = new double[scanparset_vals.size()];
double* difweights = new double[scanparset_vals.size()];
int nweights = scanparset_vals.size();
// int NParSets = scanparset_vals.size();
// Loop over all event I/O
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
for (size_t i = 0; i < eventkeys.size(); i++) {
nuiskey key = eventkeys.at(i);
// Get I/O
std::string inputfilename = key.GetS("input");
if (inputfilename.empty()) {
ERR(FTL) << "No input given for set of input events!" << std::endl;
throw;
}
std::string outputfilename = key.GetS("output");
if (outputfilename.empty()) {
outputfilename = inputfilename + ".nuisance.root";
ERR(FTL) << "No output give for set of output events! Saving to "
<< outputfilename << std::endl;
}
// Make a new input handler
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(inputfilename, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inptype =
InputUtils::ParseInputType(file_descriptor[0]);
// Make handlers for input and output
InputHandlerBase* input = InputUtils::CreateInputHandler("rawevents", inptype, file_descriptor[1]);
InputHandlerBase* output = InputUtils::CreateInputHandler("splineevents", InputUtils::kEVSPLN_Input, outputfilename);
// Get Base Events for each case.
FitEvent* rawevent = input->FirstNuisanceEvent();
FitEvent* splevent = output->FirstNuisanceEvent();
// Setup outputfile
std::string outputtest = outputfilename + ".splinetest.1DEventScan.root";
TFile* outputtestfile = new TFile(outputtest.c_str(), "RECREATE");
outputtestfile->cd();
// Save Parameter Sets
for (size_t i = 0; i < scanparset_hists.size(); i++) {
scanparset_hists[i]->Write(Form("Paramater_Set_%i", (int)i));
}
// Save a TTree of weights and differences.
TTree* weighttree = new TTree("weightscan", "weightscan");
// Make a branch for each weight set
for (size_t i = 0; i < scanparset_hists.size(); i++) {
weighttree->Branch(Form("RawWeights_Set_%i", (int)i), &rawweights[i], Form("RawWeights_Set_%i/D", (int)i) );
weighttree->Branch(Form("SplineWeights_Set_%i", (int)i), &splweights[i], Form("SplineWeights_Set_%i/D", (int)i) );
weighttree->Branch(Form("DifWeights_Set_%i", (int)i), &difweights[i], Form("DifWeights_Set_%i/D", (int)i) );
}
// Count
// int i = 0;
int nevents = input->GetNEvents();
int lasttime = time(NULL);
// Load N Chunks of the Weights into Memory
// Split into N processing chunks
int nchunks = FitPar::Config().GetParI("spline_chunks");
if (nchunks <= 0) nchunks = 1;
if (nchunks >= nevents / 2) nchunks = nevents / 2;
std::cout << "Starting NChunks " << nchunks << std::endl;
for (int ichunk = 0; ichunk < nchunks; ichunk++) {
// Skip to only do one processing chunk
// if (procchunk != -1 and procchunk != ichunk) continue;
LOG(FIT) << "On Processing Chunk " << ichunk << std::endl;
int neventsinchunk = nevents / nchunks;
int loweventinchunk = neventsinchunk * ichunk;
// int higheventinchunk = neventsinchunk * (ichunk + 1);
double** allrawweights = new double*[neventsinchunk];
double** allsplweights = new double*[neventsinchunk];
double** alldifweights = new double*[neventsinchunk];
for (int k = 0; k < neventsinchunk; k++){
allrawweights[k] = new double[nweights];
allsplweights[k] = new double[nweights];
alldifweights[k] = new double[nweights];
}
// Start Set Processing Here.
for (int iset = 0; iset < nweights; iset++) {
// Reconfigure
fRW->SetAllDials(&scanparset_vals[iset][0], scanparset_vals[iset].size());
fRW->Reconfigure();
// Reconfigure spline RW
splweight->SetAllDials(&scanparset_vals[iset][0], scanparset_vals[iset].size());
splweight->Reconfigure();
splevent->fSplineRead->SetNeedsReconfigure(true);
// Could reorder this to save the weightconts in order instead of reconfiguring per event.
// Loop over all events and fill the TTree
for (int i = 0; i < neventsinchunk; i++){
rawevent = input->GetNuisanceEvent(i+loweventinchunk);
splevent = output->GetNuisanceEvent(i+loweventinchunk);
allrawweights[i][iset] = fRW->CalcWeight(rawevent);
allsplweights[i][iset] = splweight->CalcWeight(splevent);
alldifweights[i][iset] = allsplweights[i][iset] - allrawweights[i][iset];
}
std::ostringstream timestring;
int timeelapsed = time(NULL) - lasttime;
if (timeelapsed) {
lasttime = time(NULL);
int setsleft = (nweights-iset-1) + (nweights * (nchunks - ichunk - 1));
float proj = (float(setsleft) *timeelapsed) / 60 / 60;
timestring << setsleft << " sets remaining. Last one took " << timeelapsed << ". " << proj << " hours remaining.";
}
LOG(REC) << "Processed Set " << iset << "/" << nweights <<" in chunk " << ichunk << "/" << nchunks << " " << timestring.str() << std::endl;
}
// Fill weights for this chunk into the TTree
for (int k = 0; k < neventsinchunk; k++){
for (int l = 0; l < nweights; l++){
rawweights[l] = allrawweights[k][l];
splweights[l] = allsplweights[k][l];
difweights[l] = alldifweights[k][l];
}
weighttree->Fill();
}
}
// Loop over nchunks
// Loop over parameter sets
// Set All Dials and reconfigure
// Loop over events in chunk
// Fill Chunkweightcontainers
// Once all dials are done, fill the weight tree
// Iterator to next chunk
outputtestfile->cd();
weighttree->Write();
outputtestfile->Close();
}
}
/*
// Make a high resolution spline set.
std::vector<double> nomvals = fRW->GetDialValues();
int testres = FitPar::Config().GetParI("spline_test_resolution");
std::vector< std::vector<double> > scanparset_vals;
std::vector< TH1D* > scanparset_hists;
// Loop over all params
// Add Parameters
for (size_t i = 0; i < parameterkeys.size(); i++) {
nuiskey key = parameterkeys[i];
// Get Par Name
std::string name = key.GetS("name");
if (!key.Has("low") or !key.Has("high") or !key.Has("step")) {
continue;
}
// Push Back Scan
double low = key.GetD("low");
double high = key.GetD("high");
double cur = low;
double step = key.GetD("step");
while (cur <= high) {
// Make new set
std::vector<double> newvals = nomvals;
newvals[i] = cur;
// Add to vects
scanparset_vals.push_back(newvals);
TH1D* parhist = (TH1D*)parhisttemplate->Clone();
for (size_t j = 0; j < newvals.size(); j++) {
parhist->SetBinContent(j + 1, newvals[j]);
}
scanparset_hists.push_back(parhist);
// Move to next one
cur += step;
}
}
// Print out the parameter set to test
for (int i = 0; i < scanparset_vals.size(); i++) {
std::cout << "Parset " << i;
for (int j = 0 ; j < scanparset_vals[i].size(); j++) {
std::cout << " " << scanparset_vals[i][j];
}
std::cout << std::endl;
}
// Weight holders
double* rawweights = new double[scanparset_vals.size()];
double* splweights = new double[scanparset_vals.size()];
double* difweights = new double[scanparset_vals.size()];
int NParSets = scanparset_vals.size();
// Loop over all event I/O
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
for (size_t i = 0; i < eventkeys.size(); i++) {
nuiskey key = eventkeys.at(i);
// Get I/O
std::string inputfilename = key.GetS("input");
if (inputfilename.empty()) {
ERR(FTL) << "No input given for set of input events!" << std::endl;
throw;
}
std::string outputfilename = key.GetS("output");
if (outputfilename.empty()) {
outputfilename = inputfilename + ".nuisance.root";
ERR(FTL) << "No output give for set of output events! Saving to "
<< outputfilename << std::endl;
}
// Make a new input handler
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(inputfilename, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inptype =
InputUtils::ParseInputType(file_descriptor[0]);
// Make handlers for input and output
InputHandlerBase* input = InputUtils::CreateInputHandler("rawevents", inptype, file_descriptor[1]);
InputHandlerBase* output = InputUtils::CreateInputHandler("splineevents", InputUtils::kEVSPLN_Input, outputfilename);
// Get Base Events for each case.
FitEvent* rawevent = input->FirstNuisanceEvent();
FitEvent* splevent = output->FirstNuisanceEvent();
// Setup outputfile
std::string outputtest = outputfilename + ".splinetest.1DEventScan.root";
TFile* outputtestfile = new TFile(outputtest.c_str(), "RECREATE");
outputtestfile->cd();
// Save Parameter Sets
for (size_t i = 0; i < scanparset_hists.size(); i++) {
scanparset_hists[i]->Write(Form("Paramater_Set_%i", (int)i));
}
// Save a TTree of weights and differences.
TTree* weighttree = new TTree("weightscan", "weightscan");
// Make a branch for each weight set
for (size_t i = 0; i < scanparset_hists.size(); i++) {
weighttree->Branch(Form("RawWeights_Set_%i", (int)i), &rawweights[i], Form("RawWeights_Set_%i/D", (int)i) );
weighttree->Branch(Form("SplineWeights_Set_%i", (int)i), &splweights[i], Form("SplineWeights_Set_%i/D", (int)i) );
weighttree->Branch(Form("DifWeights_Set_%i", (int)i), &difweights[i], Form("DifWeights_Set_%i/D", (int)i) );
}
// Count
int i = 0;
int nevents = input->GetNEvents();
while (rawevent and splevent) {
// Loop over 1D parameter sets.
for (size_t j = 0; j < scanparset_vals.size(); j++) {
// Reconfigure
fRW->SetAllDials(&scanparset_vals[j][0], scanparset_vals[j].size());
fRW->Reconfigure();
// Reconfigure spline RW
splweight->SetAllDials(&scanparset_vals[j][0], scanparset_vals[j].size());
splweight->Reconfigure();
splevent->fSplineRead->SetNeedsReconfigure(true);
// Calc weight for both events
rawweights[j] = fRW->CalcWeight(rawevent);
splweights[j] = splweight->CalcWeight(splevent);
difweights[j] = splweights[j] - rawweights[j];
}
if (i % 1000 == 0) {
LOG(FIT) << "Processed " << i << "/" << nevents << std::endl;
}
// Fill Array
weighttree->Fill();
// Iterate to next event.
i++;
rawevent = input->NextNuisanceEvent();
splevent = output->NextNuisanceEvent();
}
outputtestfile->cd();
weighttree->Write();
outputtestfile->Close();
}
}
*/
//*************************************
void SplineRoutines::TestSplines_NDEventThrow() {
//*************************************
// Setup RW Engine
if (fRW) delete fRW;
SetupRWEngine();
// Make a spline RW Engine too.
FitWeight* splweight = new FitWeight("splinerwaweight");
std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
TH1D* parhisttemplate = new TH1D("parhist", "parhist", parameterkeys.size(), 0.0, float(parameterkeys.size()));
// Add Parameters
for (size_t i = 0; i < parameterkeys.size(); i++) {
nuiskey key = parameterkeys[i];
std::string parname = key.GetS("name");
std::string partype = key.GetS("type");
double nom = key.GetD("nominal");
parhisttemplate->SetBinContent(i + 1, nom);
parhisttemplate->GetXaxis()->SetBinLabel(i + 1, parname.c_str());
splweight->IncludeDial( key.GetS("name"),
kSPLINEPARAMETER, nom);
splweight->SetDialValue( key.GetS("name"), key.GetD("nominal") );
}
splweight->Reconfigure();
// Make a high resolution spline set.
std::vector<double> nomvals = fRW->GetDialValues();
// int testres = FitPar::Config().GetParI("spline_test_resolution");
std::vector< std::string > scanparset_names;
std::vector< std::vector<double> > scanparset_vals;
std::vector< TH1D* > scanparset_hists;
// Loop over all params
// Add Parameters
int nthrows = FitPar::Config().GetParI("spline_test_throws");
for (int i = 0; i < nthrows; i++) {
std::vector<double> newvals = nomvals;
for (size_t j = 0; j < parameterkeys.size(); j++) {
nuiskey key = parameterkeys[j];
if (!key.Has("low") or !key.Has("high") or !key.Has("step")) {
continue;
}
// Push Back Scan
double low = key.GetD("low");
double high = key.GetD("high");
newvals[j] = gRandom->Uniform(low, high);
}
// Add to vects
scanparset_vals.push_back(newvals);
TH1D* parhist = (TH1D*)parhisttemplate->Clone();
for (size_t j = 0; j < newvals.size(); j++) {
parhist->SetBinContent(j + 1, newvals[j]);
}
scanparset_hists.push_back(parhist);
}
// Print out the parameter set to test
for (uint i = 0; i < scanparset_vals.size(); i++) {
std::cout << "Parset " << i;
for (uint j = 0 ; j < scanparset_vals[i].size(); j++) {
std::cout << " " << scanparset_vals[i][j];
}
std::cout << std::endl;
}
// Weight holders
double* rawweights = new double[scanparset_vals.size()];
double* splweights = new double[scanparset_vals.size()];
double* difweights = new double[scanparset_vals.size()];
int nweights = scanparset_vals.size();
// int NParSets = scanparset_vals.size();
// Loop over all event I/O
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
for (size_t i = 0; i < eventkeys.size(); i++) {
nuiskey key = eventkeys.at(i);
// Get I/O
std::string inputfilename = key.GetS("input");
if (inputfilename.empty()) {
ERR(FTL) << "No input given for set of input events!" << std::endl;
throw;
}
std::string outputfilename = key.GetS("output");
if (outputfilename.empty()) {
outputfilename = inputfilename + ".nuisance.root";
ERR(FTL) << "No output give for set of output events! Saving to "
<< outputfilename << std::endl;
}
// Make a new input handler
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(inputfilename, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inptype =
InputUtils::ParseInputType(file_descriptor[0]);
// Make handlers for input and output
InputHandlerBase* input = InputUtils::CreateInputHandler("rawevents", inptype, file_descriptor[1]);
InputHandlerBase* output = InputUtils::CreateInputHandler("splineevents", InputUtils::kEVSPLN_Input, outputfilename);
// Get Base Events for each case.
FitEvent* rawevent = input->FirstNuisanceEvent();
FitEvent* splevent = output->FirstNuisanceEvent();
// Setup outputfile
std::string outputtest = outputfilename + ".splinetest.NDEventThrow.root";
TFile* outputtestfile = new TFile(outputtest.c_str(), "RECREATE");
outputtestfile->cd();
// Save Parameter Sets
for (size_t i = 0; i < scanparset_hists.size(); i++) {
scanparset_hists[i]->Write(Form("Paramater_Set_%i", (int)i));
}
// Save a TTree of weights and differences.
TTree* weighttree = new TTree("weightscan", "weightscan");
// Make a branch for each weight set
for (size_t i = 0; i < scanparset_hists.size(); i++) {
weighttree->Branch(Form("RawWeights_Set_%i", (int)i), &rawweights[i], Form("RawWeights_Set_%i/D", (int)i) );
weighttree->Branch(Form("SplineWeights_Set_%i", (int)i), &splweights[i], Form("SplineWeights_Set_%i/D", (int)i) );
weighttree->Branch(Form("DifWeights_Set_%i", (int)i), &difweights[i], Form("DifWeights_Set_%i/D", (int)i) );
}
// Count
// int i = 0;
int nevents = input->GetNEvents();
int lasttime = time(NULL);
// Load N Chunks of the Weights into Memory
// Split into N processing chunks
int nchunks = FitPar::Config().GetParI("spline_chunks");
if (nchunks <= 0) nchunks = 1;
if (nchunks >= nevents / 2) nchunks = nevents / 2;
std::cout << "Starting NChunks " << nchunks << std::endl;
for (int ichunk = 0; ichunk < nchunks; ichunk++) {
// Skip to only do one processing chunk
// if (procchunk != -1 and procchunk != ichunk) continue;
LOG(FIT) << "On Processing Chunk " << ichunk << std::endl;
int neventsinchunk = nevents / nchunks;
int loweventinchunk = neventsinchunk * ichunk;
// int higheventinchunk = neventsinchunk * (ichunk + 1);
double** allrawweights = new double*[neventsinchunk];
double** allsplweights = new double*[neventsinchunk];
double** alldifweights = new double*[neventsinchunk];
for (int k = 0; k < neventsinchunk; k++){
allrawweights[k] = new double[nweights];
allsplweights[k] = new double[nweights];
alldifweights[k] = new double[nweights];
}
// Start Set Processing Here.
for (int iset = 0; iset < nweights; iset++) {
// Reconfigure
fRW->SetAllDials(&scanparset_vals[iset][0], scanparset_vals[iset].size());
fRW->Reconfigure();
// Reconfigure spline RW
splweight->SetAllDials(&scanparset_vals[iset][0], scanparset_vals[iset].size());
splweight->Reconfigure();
splevent->fSplineRead->SetNeedsReconfigure(true);
// Could reorder this to save the weightconts in order instead of reconfiguring per event.
// Loop over all events and fill the TTree
for (int i = 0; i < neventsinchunk; i++){
rawevent = input->GetNuisanceEvent(i+loweventinchunk);
splevent = output->GetNuisanceEvent(i+loweventinchunk);
allrawweights[i][iset] = fRW->CalcWeight(rawevent);
allsplweights[i][iset] = splweight->CalcWeight(splevent);
alldifweights[i][iset] = allsplweights[i][iset] - allrawweights[i][iset];
}
std::ostringstream timestring;
int timeelapsed = time(NULL) - lasttime;
if (timeelapsed) {
lasttime = time(NULL);
int setsleft = (nweights-iset-1) + (nweights * (nchunks - ichunk - 1));
float proj = (float(setsleft) *timeelapsed) / 60 / 60;
timestring << setsleft << " sets remaining. Last one took " << timeelapsed << ". " << proj << " hours remaining.";
}
LOG(REC) << "Processed Set " << iset << "/" << nweights <<" in chunk " << ichunk << "/" << nchunks << " " << timestring.str() << std::endl;
}
// Fill weights for this chunk into the TTree
for (int k = 0; k < neventsinchunk; k++){
for (int l = 0; l < nweights; l++){
rawweights[l] = allrawweights[k][l];
splweights[l] = allsplweights[k][l];
difweights[l] = alldifweights[k][l];
}
weighttree->Fill();
}
}
// Loop over nchunks
// Loop over parameter sets
// Set All Dials and reconfigure
// Loop over events in chunk
// Fill Chunkweightcontainers
// Once all dials are done, fill the weight tree
// Iterator to next chunk
outputtestfile->cd();
weighttree->Write();
outputtestfile->Close();
}
}
void SplineRoutines::SaveSplinePlots() {
if (fRW) delete fRW;
SetupRWEngine();
// Setup the spline reader
SplineWriter* splwrite = new SplineWriter(fRW);
std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
// Add splines to splinewriter
for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
iter != splinekeys.end(); iter++) {
nuiskey splkey = (*iter);
// Add Spline Info To Reader
splwrite->AddSpline(splkey);
}
splwrite->SetupSplineSet();
// Event Loop
// Loop over all events and calculate weights for each parameter set.
// Generate a set of nominal events
// Method, Loop over inputs, create input handler, then create a ttree
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
for (size_t i = 0; i < eventkeys.size(); i++) {
nuiskey key = eventkeys.at(i);
// Get I/O
std::string inputfilename = key.GetS("input");
if (inputfilename.empty()) {
ERR(FTL) << "No input given for set of input events!" << std::endl;
throw;
}
std::string outputfilename = key.GetS("output");
if (outputfilename.empty()) {
outputfilename = inputfilename + ".nuisance.root";
ERR(FTL) << "No output give for set of output events! Saving to "
<< outputfilename << std::endl;
}
// Make new outputfile
outputfilename += ".SplinePlots.root";
TFile* outputfile = new TFile(outputfilename.c_str(), "RECREATE");
outputfile->cd();
// Make a new input handler
std::vector<std::string> file_descriptor =
GeneralUtils::ParseToStr(inputfilename, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inptype =
InputUtils::ParseInputType(file_descriptor[0]);
InputHandlerBase* input = InputUtils::CreateInputHandler("eventsaver", inptype, file_descriptor[1]);
// Get info from inputhandler
int nevents = input->GetNEvents();
int countwidth = (nevents / 1000);
FitEvent* nuisevent = input->FirstNuisanceEvent();
outputfile->cd();
// int lasttime = time(NULL);
TCanvas* fitcanvas = NULL;
// Loop over all events and fill the TTree
while (nuisevent) {
// std::cout << "Fitting event " << i << std::endl;
// Calculate the weights for each parameter set
splwrite->GetWeightsForEvent(nuisevent);
splwrite->FitSplinesForEvent(fitcanvas, true);
if (fitcanvas) {
outputfile->cd();
fitcanvas->Write(Form("Event_SplineCanvas_%i", (int)i));
}
// Logging
if (i % countwidth == 0) {
LOG(REC) << "Saved " << i << "/" << nevents << " nuisance spline plots. " << std::endl;
}
// Iterate
i++;
nuisevent = input->NextNuisanceEvent();
}
// Save flux and close file
outputfile->cd();
// Close Output
outputfile->Close();
// Delete Inputs
delete input;
}
// remove Keys
eventkeys.clear();
}
void SplineRoutines::TestSplines_NDLikelihoodThrow() {
// Setup RW Engine
if (fRW) delete fRW;
SetupRWEngine();
// Make a spline RW Engine too.
FitWeight* splweight = new FitWeight("splinerwaweight");
std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
TH1D* parhisttemplate = new TH1D("parhist", "parhist", parameterkeys.size(), 0.0, float(parameterkeys.size()));
// Add Parameters
for (size_t i = 0; i < parameterkeys.size(); i++) {
nuiskey key = parameterkeys[i];
std::string parname = key.GetS("name");
std::string partype = key.GetS("type");
double nom = key.GetD("nominal");
parhisttemplate->SetBinContent(i + 1, nom);
parhisttemplate->GetXaxis()->SetBinLabel(i + 1, parname.c_str());
splweight->IncludeDial( key.GetS("name"),
kSPLINEPARAMETER, nom);
splweight->SetDialValue( key.GetS("name"), key.GetD("nominal") );
}
splweight->Reconfigure();
// Make a high resolution spline set.
std::vector<double> nomvals = fRW->GetDialValues();
// int testres = FitPar::Config().GetParI("spline_test_resolution");
std::vector< std::string > scanparset_names;
std::vector< std::vector<double> > scanparset_vals;
std::vector< TH1D* > scanparset_hists;
// Loop over all params
// Add Parameters
int nthrows = FitPar::Config().GetParI("spline_test_throws");
for (int i = 0; i < nthrows; i++) {
std::vector<double> newvals = nomvals;
for (size_t j = 0; j < parameterkeys.size(); j++) {
nuiskey key = parameterkeys[j];
if (!key.Has("low") or !key.Has("high") or !key.Has("step")) {
continue;
}
// Push Back Scan
double low = key.GetD("low");
double high = key.GetD("high");
newvals[j] = gRandom->Uniform(low, high);
}
// Add to vects
scanparset_vals.push_back(newvals);
TH1D* parhist = (TH1D*)parhisttemplate->Clone();
for (size_t j = 0; j < newvals.size(); j++) {
parhist->SetBinContent(j + 1, newvals[j]);
}
scanparset_hists.push_back(parhist);
}
// Print out the parameter set to test
for (uint i = 0; i < scanparset_vals.size(); i++) {
std::cout << "Parset " << i;
for (uint j = 0 ; j < scanparset_vals[i].size(); j++) {
std::cout << " " << scanparset_vals[i][j];
}
std::cout << std::endl;
}
// Make a new set of Raw/Spline Sample Keys
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
std::vector<nuiskey> testkeys = Config::QueryKeys("sampletest");
std::vector<nuiskey> rawkeys;
std::vector<nuiskey> splkeys;
for (std::vector<nuiskey>::iterator iter = testkeys.begin();
iter != testkeys.end(); iter++) {
nuiskey key = (*iter);
std::string samplename = key.GetS("name");
std::string eventsid = key.GetS("inputid");
nuiskey eventskey = Config::QueryLastKey("events", "id=" + eventsid);
std::string rawfile = eventskey.GetS("input");
std::string splfile = eventskey.GetS("output");
nuiskey rawkeytemp = Config::CreateKey("sample");
rawkeytemp.SetS("name", samplename);
rawkeytemp.SetS("input", rawfile);
nuiskey splkeytemp = Config::CreateKey("sample");
splkeytemp.SetS("name", samplename);
splkeytemp.SetS("input", "EVSPLN:" + splfile);
rawkeys.push_back(rawkeytemp);
splkeys.push_back(splkeytemp);
}
if (fOutputRootFile) delete fOutputRootFile;
fOutputRootFile = new TFile(fOutputFile.c_str(), "RECREATE");
fOutputRootFile->ls();
// Make two new JointFCN
JointFCN* rawfcn = new JointFCN(rawkeys, fOutputRootFile);
JointFCN* splfcn = new JointFCN(splkeys, fOutputRootFile);
// Create iteration tree in output file
fOutputRootFile->cd();
rawfcn->CreateIterationTree("raw_iterations", fRW);
splfcn->CreateIterationTree("spl_iterations", splweight);
// Loop over parameter sets.
for (size_t j = 0; j < scanparset_vals.size(); j++) {
FitBase::SetRW(fRW);
double rawtotal = rawfcn->DoEval(&scanparset_vals[j][0]);
FitBase::SetRW(splweight);
double spltotal = splfcn->DoEval(&scanparset_vals[j][0]);
LOG(FIT) << "RAW SPLINE DIF = " << rawtotal << " " << spltotal << " " << spltotal - rawtotal << std::endl;
}
fOutputRootFile->cd();
rawfcn->WriteIterationTree();
splfcn->WriteIterationTree();
}
void SplineRoutines::TestSplines_1DLikelihoodScan() {
// Setup RW Engine.
if (fRW) delete fRW;
SetupRWEngine();
// Setup Parameter Set.
// Make a spline RW Engine too.
FitWeight* splweight = new FitWeight("splinerwaweight");
// std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
TH1D* parhisttemplate = new TH1D("parhist", "parhist", parameterkeys.size(), 0.0, float(parameterkeys.size()));
// Add Parameters
for (size_t i = 0; i < parameterkeys.size(); i++) {
nuiskey key = parameterkeys[i];
std::string parname = key.GetS("name");
std::string partype = key.GetS("type");
double nom = key.GetD("nominal");
parhisttemplate->SetBinContent(i + 1, nom);
parhisttemplate->GetXaxis()->SetBinLabel(i + 1, parname.c_str());
splweight->IncludeDial( key.GetS("name"),
kSPLINEPARAMETER, nom);
splweight->SetDialValue( key.GetS("name"), key.GetD("nominal") );
}
splweight->Reconfigure();
// Make a high resolution spline set.
std::vector<double> nomvals = fRW->GetDialValues();
// int testres = FitPar::Config().GetParI("spline_test_resolution");
std::vector< std::vector<double> > scanparset_vals;
std::vector< TH1D* > scanparset_hists;
// Loop over all params
// Add Parameters
for (size_t i = 0; i < parameterkeys.size(); i++) {
nuiskey key = parameterkeys[i];
// Get Par Name
std::string name = key.GetS("name");
if (!key.Has("low") or !key.Has("high") or !key.Has("step")) {
continue;
}
// Push Back Scan
double low = key.GetD("low");
double high = key.GetD("high");
double cur = low;
double step = key.GetD("step");
while (cur <= high) {
// Make new set
std::vector<double> newvals = nomvals;
newvals[i] = cur;
// Add to vects
scanparset_vals.push_back(newvals);
TH1D* parhist = (TH1D*)parhisttemplate->Clone();
for (size_t j = 0; j < newvals.size(); j++) {
parhist->SetBinContent(j + 1, newvals[j]);
}
scanparset_hists.push_back(parhist);
// Move to next one
cur += step;
}
}
// Print out the parameter set to test
for (uint i = 0; i < scanparset_vals.size(); i++) {
std::cout << "Parset " << i;
for (uint j = 0 ; j < scanparset_vals[i].size(); j++) {
std::cout << " " << scanparset_vals[i][j];
}
std::cout << std::endl;
}
// Make a new set of Raw/Spline Sample Keys
std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
std::vector<nuiskey> testkeys = Config::QueryKeys("sampletest");
std::vector<nuiskey> rawkeys;
std::vector<nuiskey> splkeys;
for (std::vector<nuiskey>::iterator iter = testkeys.begin();
iter != testkeys.end(); iter++) {
nuiskey key = (*iter);
std::string samplename = key.GetS("name");
std::string eventsid = key.GetS("inputid");
nuiskey eventskey = Config::QueryLastKey("events", "id=" + eventsid);
std::string rawfile = eventskey.GetS("input");
std::string splfile = eventskey.GetS("output");
nuiskey rawkeytemp = Config::CreateKey("sample");
rawkeytemp.SetS("name", samplename);
rawkeytemp.SetS("input", rawfile);
nuiskey splkeytemp = Config::CreateKey("sample");
splkeytemp.SetS("name", samplename);
splkeytemp.SetS("input", "EVSPLN:" + splfile);
rawkeys.push_back(rawkeytemp);
splkeys.push_back(splkeytemp);
}
if (fOutputRootFile) delete fOutputRootFile;
fOutputRootFile = new TFile(fOutputFile.c_str(), "RECREATE");
fOutputRootFile->ls();
// Make two new JointFCN
JointFCN* rawfcn = new JointFCN(rawkeys, fOutputRootFile);
JointFCN* splfcn = new JointFCN(splkeys, fOutputRootFile);
// Create iteration tree in output file
fOutputRootFile->cd();
rawfcn->CreateIterationTree("raw_iterations", fRW);
splfcn->CreateIterationTree("spl_iterations", splweight);
// Loop over parameter sets.
for (size_t j = 0; j < scanparset_vals.size(); j++) {
FitBase::SetRW(fRW);
double rawtotal = rawfcn->DoEval(&scanparset_vals[j][0]);
FitBase::SetRW(splweight);
double spltotal = splfcn->DoEval(&scanparset_vals[j][0]);
LOG(FIT) << "RAW SPLINE DIF = " << rawtotal << " " << spltotal << " " << spltotal - rawtotal << std::endl;
}
fOutputRootFile->cd();
rawfcn->WriteIterationTree();
splfcn->WriteIterationTree();
}
/*
MISC Functions
*/
//*************************************
int SplineRoutines::GetStatus() {
//*************************************
return 0;
}