Project

General

Profile

Feature #16129

Can't find nearest wire for position

Added by Tingjun Yang over 3 years ago. Updated over 3 years ago.

Status:
Closed
Priority:
Normal
Category:
Simulation
Target version:
-
Start date:
04/06/2017
Due date:
% Done:

100%

Estimated time:
6.00 h
Spent time:
Experiment:
DUNE
Co-Assignees:
Duration:

Description

Leigh Whitehead reported the following warning:

I was running some g4 steps for protoDUNE-SP today in v06_29_00 and saw some error messages that I hadn’t seen before. Maybe they are related to Gianluca’s recent changes?

%MSG-w LArVoxelReadout:  LArG4:largeant 04-Apr-2017 08:33:45 CDT  run: 1 subRun: 1 event: 10
unable to drift electrons from point (-279.785,343.693,695.005) with exception ---- Geometry BEGIN
  Can't find nearest wire for position ( -359.416 ; 343.849 ; 694.986 ) in plane C:0 T:9 P:2 approx wire number # 479 (capped from 481)
---- Geometry END

%MSG

This needs to be understood.

Associated revisions

Revision 360f19cf (diff)
Added by Gianluca Petrillo over 3 years ago

Added option "ChargeRecoveryMargin" to LArG4.

This is related to issue #16129.
See larg4::LArVoxelReadout::RecoverOffPlaneDeposit() documentation for
details.
The projection area is the volume the plane is defined by, from the GDML
description.

Revision dab9bc6b (diff)
Added by Gianluca Petrillo over 3 years ago

Added option "ChargeRecoveryMargin" to LArG4.

This is related to issue #16129.
See larg4::LArVoxelReadout::RecoverOffPlaneDeposit() documentation for
details.
The projection area is the volume the plane is defined by, from the GDML
description.

Revision f56014d5 (diff)
Added by Gianluca Petrillo over 3 years ago

Added the concept of "active plane" in geo::PlaneGeo.

See geo::PlaneGeo::DeltaFromActivePlane() documentation for details.
This is instrumental to the solution of issue #16129.

Revision b8c6f052 (diff)
Added by Gianluca Petrillo over 3 years ago

Added the concept of "active plane" in geo::PlaneGeo.

See geo::PlaneGeo::DeltaFromActivePlane() documentation for details.
This is instrumental to the solution of issue #16129.

Revision 25c438ea (diff)
Added by Gianluca Petrillo over 3 years ago

Use the active plane as boundary for charge recovery.

This should complete resolution of issue #16129.

Revision 0afcc037 (diff)
Added by Gianluca Petrillo over 3 years ago

Use the active plane as boundary for charge recovery.

This should complete resolution of issue #16129.

Revision ad863488 (diff)
Added by Gianluca Petrillo over 3 years ago

Refined the algorithm to determine the "active area" of a wire plane.

This should provide the best solution to issue #16129 (barred bugs).

Revision a25f7ec2 (diff)
Added by Gianluca Petrillo over 3 years ago

Refined the algorithm to determine the "active area" of a wire plane.

This should provide the best solution to issue #16129 (barred bugs).

Revision 1bf8a304 (diff)
Added by Gianluca Petrillo over 3 years ago

Adding a small margin to the charge recovery position.

This is supposed to work around rounding errors.
Context is issue #16129.

Revision cb8c2a08 (diff)
Added by Gianluca Petrillo over 3 years ago

Adding a small margin to the charge recovery position.

This is supposed to work around rounding errors.
Context is issue #16129.

History

#1 Updated by Gianluca Petrillo over 3 years ago

It may be.
I don't have the geometry information in front of me, but let's say that the coverage of the wire plane is $z \in \left[ 695,995 \
ight] \textnormal{cm}$.
The charge deposited at z = 695.005 gets diffused. Some of it gets to z = 694.986 cm. That is outside of the wire plane, hence the error.

Does that make sense?

If the old ChannelMappingAlg was capping to the "closest" wire, you would not have noticed that.

#2 Updated by Thomas Junk over 3 years ago

Yes, in DuneApaChannelMapAlg.cxx, in the method:

DuneApaChannelMapAlg::NearestWireID

the coordinates are silently adjusted by pushing Y and Z (gasp, hard-coded Y and Z!) so that
they are inside the range covered by the wires. The assumption here is that drifting charge
diverts inwards towards the collection wires along Z if it is outside the Z range, and along Y if it is outside the Y range.

But this method is no longer called (I put breakpoints on it). Instead, in PlaneGeo.cxx,

PlaneGeo::NearestWireID

caps the wire ID and does so noisily by throwing an exception. If you look at the error messages, they mostly are small nudges, by one wire or two. I assume that if the exception is thrown, then the charge is not drifted so we lose some calorimetry info here and really should fix it.

I see some methods that look promising in PlaneGeo.cxx -- PlaneGeo::DeltaFromPlane does some calling of symmetricCapDelta. Are these caps being set right? In the DUNE caps, I looked at the wire endpoints. If you use the frame size, you might get some space off the ends as the wires don't got all the way out to the ends of the frames? Especially in DUNE, the wires wrap around FR-4 boards with different positions so it depends on which plane you're in.

#3 Updated by Thomas Junk over 3 years ago

Our preference would be to keep our old behavior -- cap Y and Z so that drifting charge lands in the part of the plane covered by the wires. Silently capping the wire number and not throwing an exception is an improvement but still problematic as the diversion of the charge goes in the wrong direction. This is a bigger problem if the active volume is much bigger than the wire planes, but we have had a presentation where a user noticed the difference even with the geometry as is.

#4 Updated by Gianluca Petrillo over 3 years ago

  • Assignee changed from Tingjun Yang to Gianluca Petrillo
  • Estimated time set to 6.00 h

I am going to steal this from Tingjun.

I believe the methods Tom mentions use the dimensions of the volume which defines the wire plane, volPlaneXxx.
I don't know if that includes the wire plane frame.

Since the methods in geo::PlaneGeo are designed for this, you get to decide how they should behave.

One option is to trust the detector description and use volPlaneXxx directly. The other option is to use the coverage dimensions.
I think you mean you want the latter, but I would like a confirmation.
Since this will likely noticeably slow down the simulation, it will be optional and disabled by default.
Task for you: find a good name for the configuration parameter controlling this feature.

#5 Updated by Thomas Junk over 3 years ago

The way we handled this in dunetpc was on initialization, for each TPC, the endpoints of all wires were queried, and the max and min Y and Z saved -- see the Initialize method of DuneApaChannelMapAlg.cxx. I guess this would have to generalize in larsoft if we're not using our method anymore.

#6 Updated by Gianluca Petrillo over 3 years ago

Can you provide me with suitable input to reproduce the issue, and a description of what you need to consider the issue solved?

My suggestion for these are:
  • point me to a generated input file, and to a FHiCL configuration running LArG4
  • accept as an indication of the solution the disappearance of those warnings

Particularly important is the first item, with an event known to show the issue. The simplest the event, the better (e.g. single muon?).

#7 Updated by Thomas Junk over 3 years ago

How to reproduce:

I ran a short MARLEY sample which made some of these warnings and put the gen file (it's small)
on the app disk. The second event produces the warnings.

To reproduce from a bare login:

$ source /grid/fermiapp/products/dune/setup_dune.sh
$ setup dunetpc v06_31_00 -q debug:e10
$ lar --nskip 1 -n 1 -c standard_g4_dune10kt_1x2x6.fcl /dune/app/users/trj/jtest/gen.root

#8 Updated by Gianluca Petrillo over 3 years ago

  • Project changed from dunetpc to LArSoft
  • Category set to Simulation
  • Status changed from Assigned to Work in progress
  • % Done changed from 0 to 10
  • Experiment DUNE added

Request moved to LArSoft tracker since it involves generic code.

#9 Updated by Gianluca Petrillo over 3 years ago

  • Status changed from Work in progress to Feedback
  • % Done changed from 10 to 100

I have prepared the code which should act as the desired.
The charge is moved on the "active" part of the plane if its landing position on the plane is less than a margin from its limits.
The active area is defined as the smallest rectangle containing all the wire ends, plus half a wire pitch if the wires are aligned to a the frame side. The active area is defined by the plane object (geo::PlaneGeo).

The code is in feature/gp_Issue16129 branch of larcore and larsim. Both are necessary for the latest larsim to work.
This changes need physics testing and validation.
The feature can be enabled by setting a margin larger than 0, for example:

physics.producers.largeant.ChargeRecoveryMargin: 1.0 # cm
If the parameter is left untouched to the default of 0, the new code is completely skipped (except that at the beginning of the job plane objects always compute their active area). There is currently a single margin value, shared by all planes, all TPCs and both "width" and "depth" directions.

I could not notice any slowdown of the code when turning on the new feature, although there is certainly some since it's computing more.

I wait the feedback from your testing before declaring this resolved and ready to go into the release.

#10 Updated by Leigh Whitehead over 3 years ago

I just tried to test this with protoDUNE-SP, but I still see the error messages. I set the margin value to 1cm, 10cm, and even 1m. I assume that I shouldn't be seeing the error messages any more?

I have checked out both of your feature branches in larcore and larsim, and the develop of dunetpc in order to change the .fcl file to include the margin parameter.

#11 Updated by Leigh Whitehead over 3 years ago

Just to clarify, if I raise ChargeRecoveryMargin from 1.0cm to 2.5cm, I don't have any instances failing this block of code:

// won't recover: too far
if ((std::abs(offPlane.X()) > fOffPlaneMargin)
   || (std::abs(offPlane.Y()) > fOffPlaneMargin))
   return pos;

I still see the error messages, however.

#12 Updated by Gianluca Petrillo over 3 years ago

  • Status changed from Feedback to Assigned

Can you show me how to reproduce that, with a proper input file?

It might be that the "active" area as determined by geo::PlaneGeo is too large. Then, the points are brought to the border of that area, but that border is still too far and triggers the exception.
But I need to see it running.

#13 Updated by Leigh Whitehead over 3 years ago

Here is the input file: /pnfs/dune/scratch/users/lwhite86/v06_29_00/gen/mcc9Test/18612787_0/gen_protoDune_beam_p3GeV_cosmics_aa48f4c5-1d33-44dc-b0b1-5ddbccaa4a75.root

I used protoDUNE_g4_3ms.fcl, which calls protoDUNE_g4.fcl, to which I added the line physics.producers.largeant.ChargeRecoveryMargin: 2.5 # cm

#14 Updated by Gianluca Petrillo over 3 years ago

  • Status changed from Assigned to Feedback

I have identified the problem in rounding errors. I am now using an hard-coded margin to move the points slightly within the plane rather than exactly at the border (which turned out to be not that exactly after all).
I have also changed the code not to be subject to that rounding error any more... or so I think.

My current solution is not perfect. Charge is still lost in a small region of the planes. More precisely, it is lost around the two corners close to the two shortest wires (it is expected not to be a problem for planes with wires parallel to the frame).
In those corners, the charge half a pitch away from the last wire can still be lost (the margin value does not affect this fact).

If this is not acceptable, I will need some more invasive modifications to the code.
Please advise.

#15 Updated by Leigh Whitehead over 3 years ago

Gianluca, can you comment where exactly the charge was being lost from before your latest change? Roughly how much of a reduction in the problem are we talking about now?

#16 Updated by Gianluca Petrillo over 3 years ago

  • % Done changed from 100 to 90

Behold my brand-new wire plane:

      |     |      |      |      |      |   xxxxxxxxxx
      |     |      |      |      |      |   xxxxxxxxxx
      |     v      v      v      v      v   xxxxxxxxxx
 - - -+=====================================xxxxxxxxxx
      |`-._    `-._    `-._    `-._    `-._   xxxxxxxx
~~~~~>|    `-._    `-._    `-._    `-._    `-._|
      |`-._    `-._    `-._    `-._    `-._    |<~~~~~
~~~~~>|    `-._    `-._    `-._    `-._    `-._|
      |`-._    `-._    `-._    `-._    `-._    |<~~~~~
~~~~~>|    `-._    `-._    `-._    `-._    `-._|
      |`-._    `-._    `-._    `-._    `-._    |<~~~~~
xxxxxxxx   `-._    `-._    `-._    `-._    `-._|
xxxxxxxxxx=====================================+- - - 
xxxxxxxxxx   ^      ^      ^      ^      ^     |
xxxxxxxxxx   |      |      |      |      |     |
xxxxxxxxxx   |      |      |      |      |     |

The areas with x is where charge is lost. That area starts half a wire pitch away from the last wire.
Which is in fact pretty much as it was with the exceptions, in the end. The difference is in the side areas, which are now treated kind-of-correctly (where the arrows are, which remind us that the charge is drifted back to the border of the plane).
There should be no loss on planes with vertical wires... I think.

#17 Updated by Gianluca Petrillo over 3 years ago

  • Status changed from Feedback to Resolved
  • % Done changed from 90 to 100

I have refined the algorithm.
The idea is that there is an "active" area of the plane, and charge which is outside that area but within the specified margin is shifted to the border of the active area.
Then the simulation proceeds normally (and note that this happens after any displacement from space charge).

The problem is how to define the active area. If it is too small, physics will be somehow affected since the ends of the wires will be ignored. If it is too large, the border of the active area may be far enough from the wires, that the exception will still be triggered (which is in fact what Leigh observed).

So I spent good part of yesterday and today to develop an algorithm to determine the largest active area which includes no point farther than half a wire pitch from all wires. The tricky part is to have it for all orientation of planes, and all orientation of wires.
In the end, it seems to work and I don't see exceptions thrown.
It is pushed in the branches feature/gp_Issue16129 of larcore and larsim, as before.
Please validate its physics...

#18 Updated by Gianluca Petrillo over 3 years ago

A barely related point: for your production, I recommend against using the random policy for seeds of NuRandomService.
The perEvent policy should provide enough randomness and makes reproducing errors less of a pain.

#19 Updated by Tingjun Yang over 3 years ago

Hi Gianluca,

Could you merge develop into both branches? There seem to be conflicts.

Thanks,
Tingjun
Gianluca Petrillo wrote:

I have refined the algorithm.
The idea is that there is an "active" area of the plane, and charge which is outside that area but within the specified margin is shifted to the border of the active area.
Then the simulation proceeds normally (and note that this happens after any displacement from space charge).

The problem is how to define the active area. If it is too small, physics will be somehow affected since the ends of the wires will be ignored. If it is too large, the border of the active area may be far enough from the wires, that the exception will still be triggered (which is in fact what Leigh observed).

So I spent good part of yesterday and today to develop an algorithm to determine the largest active area which includes no point farther than half a wire pitch from all wires. The tricky part is to have it for all orientation of planes, and all orientation of wires.
In the end, it seems to work and I don't see exceptions thrown.
It is pushed in the branches feature/gp_Issue16129 of larcore and larsim, as before.
Please validate its physics...

#20 Updated by Tingjun Yang over 3 years ago

On a related issue, I now see the similar exception throw in the standard dune reconstruction chain:

%MSG-e pma::Track3D:  PMAlgTrackMaker:pmtracktc 15-Apr-2017 11:16:04 CDT  run: 20000002 subRun: 1 event: 8                                                                      
Flip, endpoint closer to vStart.                                                                                                                                                
%MSG                                                                                                                                                                            
%MSG-s ArtException:  PostProcessPath reco 15-Apr-2017 11:16:12 CDT  PostProcessEvent                                                                                           
cet::exception caught in art                                                                                                                                                    
---- EventProcessorFailure BEGIN                                                                                                                                                
  An exception occurred during current event processing
  ---- ScheduleExecutionFailure BEGIN
    ProcessingStopped.

    ---- Geometry BEGIN
      Can't find nearest wire for position ( -75.3758 ; -542.249 ; 1393.52 ) in plane C:0 T:20 P:2 approx wire number # 479 (capped from 482)
      cet::exception going through module
    ---- Geometry END
    Exception going through path reco
  ---- ScheduleExecutionFailure END
---- EventProcessorFailure END
%MSG
%MSG-s ArtException:  PostProcessPath reco 15-Apr-2017 11:16:12 CDT  PostProcessEvent
cet::exception caught in art
---- EventProcessorFailure BEGIN
  An exception occurred during current event processing
  ---- ScheduleExecutionFailure BEGIN
    ProcessingStopped.

    ---- Geometry BEGIN
      Can't find nearest wire for position ( -75.3758 ; -542.249 ; 1393.52 ) in plane C:0 T:20 P:2 approx wire number # 479 (capped from 482)
      cet::exception going through module
    ---- Geometry END
    Exception going through path reco
  ---- ScheduleExecutionFailure END
---- EventProcessorFailure END
%MSG
Art has completed and will exit with status 18.

I already checked out feature/gp_Issue16129 of larcore. Here is the command to reproduce the problem:
lar -c standard_reco_dune10kt_1x2x6.fcl /pnfs/dune/persistent/dunepro/v06_17_00/detsim/prodgenie_nue_dune10kt_1x2x6/12485578_0/prodgenie_nue_dune10kt_1x2x6_0_20161210T055906_detsim2_b2ff2530-418c-4a52-857e-44428b5e99b2.root --nskip 7 -n 1

#21 Updated by Robert Sulej over 3 years ago

Is this exception thrown when the intersection of wires in different planes cannot be found (is outside of active area of plane)?

Or it can happen when 2D projection of 3D point is calculated when it falls out of active area?

#22 Updated by Tingjun Yang over 3 years ago

Apologies. I checked out Gianluca's feature branch in larcore and built it but I did not set it up. After I set it up properly, I don't have problem with reconstruction. Sorry about the confusion.

The feature branches should be merged soon.

Tingjun
Tingjun Yang wrote:

On a related issue, I now see the similar exception throw in the standard dune reconstruction chain:
[...]

I already checked out feature/gp_Issue16129 of larcore. Here is the command to reproduce the problem:
lar -c standard_reco_dune10kt_1x2x6.fcl /pnfs/dune/persistent/dunepro/v06_17_00/detsim/prodgenie_nue_dune10kt_1x2x6/12485578_0/prodgenie_nue_dune10kt_1x2x6_0_20161210T055906_detsim2_b2ff2530-418c-4a52-857e-44428b5e99b2.root --nskip 7 -n 1

#23 Updated by Tingjun Yang over 3 years ago

  • Status changed from Resolved to Work in progress

Apologies again!

This is actually a real problem. It is in EMShowerAlg:

TVector2 shower::EMShowerAlg::Project3DPointOntoPlane(TVector3 const& point, int plane, int cryostat) {

  TVector2 wireTickPos = TVector2(-999., -999.);

  double pointPosition[3] = {point.X(), point.Y(), point.Z()};

  geo::TPCID tpcID = fGeom->FindTPCAtPosition(pointPosition);
  int tpc = 0;
  if (tpcID.isValid)
    tpc = tpcID.TPC;
  else
    return wireTickPos;

  // Construct wire ID for this point projected onto the plane
  geo::PlaneID planeID = geo::PlaneID(cryostat, tpc, plane);
  geo::WireID wireID = fGeom->NearestWireID(point, planeID);

  wireTickPos = TVector2(GlobalWire(wireID),
                         fDetProp->ConvertXToTicks(point.X(), planeID));

  // wireTickPos = TVector2(fGeom->WireCoordinate(point.Y(), point.Z(), planeID.Plane, tpc % 2, planeID.Cryostat),
  //              fDetProp->ConvertXToTicks(point.X(), planeID.Plane, tpc % 2, planeID.Cryostat));

  //return wireTickPos;
  return HitPosition(wireTickPos, planeID);

}

I will contact Mike.

Tingjun Yang wrote:

Apologies. I checked out Gianluca's feature branch in larcore and built it but I did not set it up. After I set it up properly, I don't have problem with reconstruction. Sorry about the confusion.

The feature branches should be merged soon.

Tingjun
Tingjun Yang wrote:

On a related issue, I now see the similar exception throw in the standard dune reconstruction chain:
[...]

I already checked out feature/gp_Issue16129 of larcore. Here is the command to reproduce the problem:
lar -c standard_reco_dune10kt_1x2x6.fcl /pnfs/dune/persistent/dunepro/v06_17_00/detsim/prodgenie_nue_dune10kt_1x2x6/12485578_0/prodgenie_nue_dune10kt_1x2x6_0_20161210T055906_detsim2_b2ff2530-418c-4a52-857e-44428b5e99b2.root --nskip 7 -n 1

#24 Updated by Gianluca Petrillo over 3 years ago

I have merged the current develop branches into feature/gp_Issue16129 branches of larcore and larsim.

Tingjun Yang wrote:

Apologies again!

This is actually a real problem. It is in EMShowerAlg:
[...]
I will contact Mike.

Also take a look to geo::PlaneGeo::WireCoordinate() and geo::PlaneGeo::DistanceFromPlane() methods.
These are purely projection functions which do not care about physical limits of the detector.

#25 Updated by Gianluca Petrillo over 3 years ago

  • Status changed from Work in progress to Resolved

As Leigh has approved the fix, I am marking this ticket as resolved.
The issue reported by Tingjun in note #23 is possibly related to this one, but it is not the same issue: this one is about simulation code, and the solution is in the simulation code.
If support is needed for EMShowerAlg, please contact me or open a specific tracking ticket.

#26 Updated by Tingjun Yang over 3 years ago

I used code suggested by Gianluca to patch EMShowerAlg and ClusterCrawlerAlg:
larreco:dbbeb8109a96ca923654f25b91baf28e7a992f8e and larreco:562f757d738abc4e4eee1e7525167132b1a11422

#27 Updated by Katherine Lato over 3 years ago

  • Status changed from Resolved to Closed


Also available in: Atom PDF