feat: add oscprob

models/dayabay_v0.py

@@ -20,9 +20,9 @@ def model_dayabay_v0():
                 ('i', 'isotope'): ('U235', 'U238', 'Pu239', 'Pu241'),
                 ('b', 'background'): ('acc', 'lihe', 'fastn', 'amc', 'alphan'),
                 })
-    idx_r= index.sub('r')
-    idx_rd= index.sub(('r', 'd'))
-    idx_ri= index.sub(('r', 'i'))
+    idx_r = index.sub('r')
+    idx_rd = index.sub(('r', 'd'))
+    idx_ri = index.sub(('r', 'i'))
     list_reactors = idx_r.values
     list_dr = idx_rd.values
     list_reactors_isotopes = idx_ri.values
@@ -37,6 +37,7 @@ def model_dayabay_v0():
         #
         load_parameters({'path': 'oscprob'    , 'load': datasource/'parameters/oscprob.yaml'})
         load_parameters({'path': 'oscprob'    , 'load': datasource/'parameters/oscprob_solar.yaml'})
+        load_parameters({'path': 'oscprob'    , 'load': datasource/'parameters/oscprob_constants.yaml'})
 
         load_parameters({'path': 'ibd'        , 'load': datasource/'parameters/pdg2012.yaml'})
         load_parameters({'path': 'ibd.csc'    , 'load': datasource/'parameters/ibd_constants.yaml'})
@@ -64,6 +65,7 @@ def model_dayabay_v0():
         # Create nodes
         #
         labels = LoadYaml(datasource/'labels.yaml')
+        parameters = storage('parameter')
         nodes = storage.child('nodes')
         inputs = storage.child('inputs')
         outputs = storage.child('outputs')
@@ -71,31 +73,38 @@ def model_dayabay_v0():
         from dagflow.lib.Array import Array
         from dagflow.lib.View import View
         from numpy import linspace
-        edges_costheta=Array.make_stored("edges.costheta", [-1, 1], label_from=labels)
-        edges_energy_common=Array.make_stored("edges.energy_common", linspace(0, 12, 241), label_from=labels)
-        edges_energy_enu=View.make_stored("edges.energy_enu", edges_energy_common, label_from=labels)
-        edges_energy_edep=View.make_stored("edges.energy_edep", edges_energy_common, label_from=labels)
-        edges_energy_evis=View.make_stored("edges.energy_evis", edges_energy_common, label_from=labels)
-        edges_energy_erec=View.make_stored("edges.energy_erec", edges_energy_common, label_from=labels)
-
+        edges_costheta, _=Array.make_stored("edges.costheta", [-1, 1], label_from=labels)
+        edges_energy_common, _=Array.make_stored("edges.energy_common", linspace(0, 12, 241), label_from=labels)
+        View.make_stored("edges.energy_enu", edges_energy_common, label_from=labels)
+        edges_energy_edep, _=View.make_stored("edges.energy_edep", edges_energy_common, label_from=labels)
+        View.make_stored("edges.energy_evis", edges_energy_common, label_from=labels)
+        View.make_stored("edges.energy_erec", edges_energy_common, label_from=labels)
+
+        integration_orders_edep, _=Array.from_value("integration.ordersx", 4, edges=edges_energy_edep, label_from=labels)
+        integration_orders_costheta, _=Array.from_value("integration.ordersy", 4, edges=edges_costheta, label_from=labels)
         from dagflow.lib.IntegratorGroup import IntegratorGroup
-        integration_orders_edep=Array.from_value("integration.ordersx", 4, edges=edges_energy_edep, label_from=labels)
-        integration_orders_costheta=Array.from_value("integration.ordersy", 4, edges=edges_costheta, label_from=labels)
-        integrator=IntegratorGroup.replicate('2d', 'kinematics_sampler', 'kinematics_integrator', replicate=list_dr)
+        integrator, _=IntegratorGroup.replicate('2d', 'kinematics_sampler', 'kinematics_integrator', replicate=list_dr)
         integration_orders_edep >> integrator.inputs["ordersX"]
         integration_orders_costheta >> integrator.inputs["ordersY"]
         outputs['integration.mesh_edep'] = (int_mesh_edep:=integrator.outputs['x'])
         outputs['integration.mesh_costheta'] = (int_mesh_costheta:=integrator.outputs['y'])
 
         from reactornueosc.IBDXsecO1Group import IBDXsecO1Group
-        nodes['ibd'] = (ibd:=IBDXsecO1Group(use_edep=True))
+        ibd, _ = IBDXsecO1Group.make_stored(use_edep=True)
         ibd << storage('parameter.constant.ibd')
         ibd << storage('parameter.constant.ibd.csc')
         int_mesh_edep >> ibd.inputs['edep']
         int_mesh_costheta >> ibd.inputs['costheta']
         outputs['ibd'] = ibd.outputs['result']
 
-        integrator.print()
+        from reactornueosc.NueSurvivalProbability import NueSurvivalProbability
+        NueSurvivalProbability.replicate('oscprob', replicate=list_dr)
+        ibd.outputs['enu'] >> inputs('oscprob.enu')
+        parameters('constant.baseline') >> inputs('oscprob.L')
+        nodes('oscprob') << parameters('free.oscprob')
+        nodes('oscprob') << parameters('constrained.oscprob')
+        nodes('oscprob') << parameters('constant.oscprob')
+
         ibd.outputs['result'] >> inputs('kinematics_integrator')
 
     storage('outputs').read_labels(labels, strict=True)
@@ -107,7 +116,7 @@ def model_dayabay_v0():
         print(storage.to_table(truncate=True))
         return
 
-    storage('outputs').plot(show_all=True)
+    # storage('outputs').plot(show_all=True)
 
     storage['parameter.normalized.detector.eres.b_stat'].value = 1
     storage['parameter.normalized.detector.eres.a_nonuniform'].value = 2