Plane ID in OmnibusSigProc and ROI Formation/Refinement for PDHD APA1

[wenqiang-gu]

In ProtoDUNE HD APA1, the U, W planes are basically induction, and V plane is collection. To optimize the signal processing for V & W planes, we'd like to consider them as collection and induction, respectively.

The current logic in ROI_formation and ROI_refinement rely on an input plane ID to separate induction (plane=0,1) and collection (plane=2). For example, ROI_formation::find_ROI_loose(). Similar logics are almost everywhere in ROI_formation or ROI_refinement, therefore, it appears that we need a refactorization.

To avoid a complicated refactorization, I propose a light-weight solution for PDHD APA1, which only requires: a hack in the wire geometry json, a hack in the field response (FR) json, and a minimal update in OmnibusSigProc.

1. In the wire geometry JSON, switch the order of plane ident numbers in the face declaration (e.g., 0,1,2 -> 0,2,1):

   "faces": [
     {
       "Face": {
         "ident": 0,
         "planes": [
           0,
           2,
           1
         ]
       }
     },
     {
       "Face": {
         "ident": 1,
         "planes": [
           3,
           5,
           4
         ]
       }
     },

   Such change could be passed to the construction of AnodeFace. As a result, the corresponding plane ID for U, V and W will be: OspChan::plane	IWirePlane::planeid.index() - the order in the JSON	View	IWirePlane::ident() - imported from JSON
   0	0	U	0
   1	1	W	2
   2	2	V	1

Here are more debug info printed from OmnibusSigproc when building OspChan: (plane index, anode, face, plane ident, ch ident range) Now knowing that the imported raw data m_r_data[plane] and int plane parameters in all ROI functions (e.g., find_ROI_loose(int plane, ...)) are synchronized with OspChan::plane, it would exactly meet our purpose: plane=1 for W as induction, plane=2 for V as collection.

2. In the FR JSON, hack the planeid to 2 for V, and similarly, hack it to 1 for W. Note that, this planeid is independent of other plane IDs for AnodePlane.

3. To use the hacked field response JSON in the OmnibusSigProc, we need to use fravg.plane(iplane) instead of fravg.planes[iplane]. For normal APAs, fravg.planes[iplane] and fravg.plane(iplane) returns same results. For hacked APA1 FR JSON, the fravg.plane(1) would retrutn W plane field response.

Note that, we don't need to update anything for simulation as the PlaneImpactResponse is using the fr.plane() function in the implementation.

[brettviren]

I think this "hack" is reasonable, @wenqiang-gu

[wenqiang-gu]

Thank you for the confirmation. Let me work on the PR and have some essential validation before I commit it.

[wenqiang-gu]

Here is the sigproc result with wirecell v0_27_1 (default version from dunesw v09_91_03d00) for a PDHD beam trigger:

Here is the sigproc result with wirecell v0_27_1 and the proposed hacks:

Two results are very similar, which confirms my proposal works for the signal processing. Note that, I don't expect identical results because the thresholds and filters are different.

I also wanted to test the proposal for simulation, while I realized that the WirePlaneId::index() is widely used, for example, in the Digitizer.cxx:

Accordingly, the simulated baseline would be swapped in APA1 V & W because W's wpid.index() would be 1 now.

Given that WirePlaneId::index() can be used in other places, I suggest we limit my proposal to OmnibusSigProc only. What do you think? @brettviren

[wenqiang-gu]

Just to clarify, in ProtoDUNE HD, the simulation stage and the reco stage are separate, so it's okay to configure different wire geo JSON and field response JSON for each stage.

[brettviren]

I think this looks good for fixing det+SP for the problematic APA1.

Just a comment: even in the same WCT job we may use different wires and FR for different components. Applying this "hack" does not require the "hacked" components to be isolated in their own job. Though, that isolation may be a natural consequence of other reasons.

And, limiting this "hack" to det+SP makes sense. In principle, it could be applied to sim but since sim+SP "cancel" the FR to first order it means a "hacked' sim+SP should give essentially the same results as the nominal sim+SP.

[brettviren]

@wenqiang-gu the PR reinterprets iplane from being a plane index to being a plane ident. Sorry I didn't realize this is what your item 3 above implied.

This will only work as intended "by accident" assuming all experiments define plane ident in [0,1,2]. But, experiments are free to chose numbers outside this range or even assign (nominal) U/V/W in different order than 0/1/2

Now, I believe all current experiments do in fact follow this more restrictive convention so I think this PR will work for them. But, I worry about what confusion may arise in the future.

So, I think we need to step back as we are digging a deeper hole by fixing a problem by adding another problem.

How about we do something a bit more robust? Let's ditch the "V-W swap hack" and add a new configuration item to OmnibusSigProc that explicitly tells what we want. Consider a small array named like m_is_collection that maps plane index to a boolean. The default value would be [false, false, true] while the APA1 case would require configuring [false,true,false] and we would provide the standard APA wires and FR files. Then an existing if (plane == 2) tests would become if (m_is_collection[plane]) and etc for non-collection planes.

I know this requires a bigger code change but maybe not so bad?

[wenqiang-gu]

@brettviren I think the m_is_collection trick would work for OmnibusSigproc itself, however, many functions in ROI_formation and ROI_refinement still use a plane number to specify U/V/W. I need to think about this.

[wenqiang-gu]

Hi @brettviren , after a second thought, I think it would be very difficult to swap V & W with just the m_is_collection trick. Here is a function call inside OmnibusSigProc:

roi_form.find_ROI_by_decon_itself(iplane, m_r_data[iplane], r_data_tight);

Inside ROI formation, it follows a convention: U/plane==0 , V/plane==1, W/plane==2. Perhaps we could bypass this convention using a configurable mapping inside ROI formation. However, more problem could arise because the input data m_r_data[iplane] still follow the connection from OmnibusSigProc.

Indeed, we need to avoid digging a deeper hole for other experiments. How about I add a protection for loading the field reponse? For example,

bool m_load_fr_with_plane_ident; // default false, configurable 

auto arr = Response::as_array(fravg.planes[iplane], fine_nwires, fine_nticks);
if (m_load_fr_with_plane_ident) {
    // Note, set plane ident properly in your JSON
    arr = Response::as_array(*(fravg.plane(iplane)), fine_nwires, fine_nticks);
}

[brettviren]

@wenqiang-gu I think your m_load_fr_with_plane_ident idea is a good practical solution.

Please add some code comments to describe this new config option and include a reference to this Issue.

(And, use a full if/else block in order to avoid two calls to as_array() )

Thanks for all your careful checking.

[wenqiang-gu]

Hi @brettviren , I updated my PR: https://github.com/WireCell/wire-cell-toolkit/pull/324

Could you take a look?

[wenqiang-gu]

Hi @brettviren , regarding the closed PR in larwirecell, I would like to propose a different solution for swapping V&W in OmnibusSigProc.

First of all, the drawback of the previous solution (swapping V&W in JSON) is that it changes V&W's plane.index() globally. Therefore, when larwirecell::FrameSaver saves a frame, it still uses the swapped V&W plane index: V- kWlayer, W- kVlayer.

Now I would like to propose a local hack that limits to OmnibusSigproc. In the other word, we don't need to change the JSON inputs of wire geometry or field response. Within OmnibusSigProc, there are three places that explicitly uses plane ID: OspChan::plane, m_r_data[plane], overall_resp[plane].

If we can create a configurable map between plane and the "virtual layer", we can simply change their order. For example,

Meanwhile, since the wire geometry JSON is not hacked, the function call WirePlaneId::index() would still return kVlayer for a channel number from the actual V plane, which fills the gap for the closed PR in larwirecell.

[wenqiang-gu]

Added the new solution in the PR: https://github.com/WireCell/wire-cell-toolkit/pull/325

[brettviren]

@wenqiang-gu Nice. Yes, this looks like a much more clean "hack".