Uniform random picking of points on detector surfaces for noise studies

[wdconinc]

Is there a functionality in DD4hep to pick uniformly distributed random points on a detector surface? We have one of our collaborators working on implementing random tracking pixels firing (shot noise) in addition to the simulated physics events, to study impact on track reconstruction efficiency.

I was looking at the DDRec::Surface class but don't see such functionality, but this may be implemented elsewhere. I was also looking at the DDDigi subtree (including DigiPoissonNoise) but I was not able to discern how to adapt that to this use case.

[MarkusFrankATcernch]

No there is not.
This problem is highly linked to #1297
and #580 .
A even only half way efficient generic solution is close to impossible to implement.
The impossibility of such an implementation effectively stopped generic noise handling in DDDigi.

Possible solutions are only possible if a segmentation allows to enumerate easily it's cells like rectangular wafers etc.
If you can come up with a sound idea, I am happy to implement!

[wdconinc]

Would a non-generic solution be acceptable, e.g. only certain surfaces or certain shapes at first, but with a return value that makes it clear when the request failed?

[ShujieL]

@MarkusFrankATcernch thanks for pointing us to the previous discussions.

Knowing that we can not pick random positions on surface, I tried to select random noisy pixels by cell ID. For example, for our three-layer Si vertex tracker with segmentation below:
https://github.com/eic/epic/blob/4bb200ff4e3e6f6e1b43bba55a2dc72dc9b44a4c/compact/tracking/vertex_barrel.xml#L151C1-L156C14
<id>system:8,layer:4,module:12,sensor:2,x:32:-16,y:-16</id>
I randomly generate a 64-bit cell ID, fix the last 8 bits by a given system id, then call CellIDPositionConverter to check if the cell id is valid, and accumulate number of valid noise hits:

 while (i < num_noise_hits)
  {
    uint64_t randomNum = distr(eng);                    // Generate a random 64-bit cell id
    randomNum = (randomNum & ~uint64_t(0xFF)) | vol_id; // Clear the last 8 bits and set them to vol ID(system id)
    try
    {
      auto pos = m_converter.position(randomNum);
    std::cout << std::bitset<8>(vol_id) << ":" << i << "/" << nn << ":::" << randomNum << "\n";
      std::cout << "position r,z=" << sqrt(pos.x()*pos.x() + pos.y()*pos.y()) << "  " << pos.z() << std::endl;
      i++;
    }
    catch (std::exception &e)
    {
      nn++;
    }
  }

This code does give me a list of valid cell id (despite of many volume manager complaints on invalid segmentation), and their layer and module id etc seem to be valid. However, the corresponding global position of those generated noise hits don't always make sense. For example, our barrel has radii of 3.6, 4.8, and 12cm, but the cell ID position from the converter is 30.88cm.

Any suggestions ?

[MarkusFrankATcernch]

@ShujieL For CartesianGridXY on a Box (such as a Wafer) it is actually possible to do this, because the CellIDs are easily enumerable.
How to do this, see DDDigi/plugins/CellScanner_CartesianGridXY.cpp:

    template <>
    void init_segmentation_data<CartesianGridXY>(segmentation_data<CartesianGridXY>& data,
                                                 const Segmentation& seg)

and

    template <typename SEGMENTATION, typename SOLID> void
    CellScanner<SEGMENTATION,SOLID>::operator()(DigiContext& context,
                                                PlacedVolume pv,
                                                VolumeID vid,
                                                const cell_handler_t& cell_handler)

at line 58:

      long nx = 2e0 * b->GetDX() / segment.x_grid_size;
      long ny = 2e0 * b->GetDY() / segment.y_grid_size;

instead of looping over all possible cell ids like in this example,
simply throw a dice and then you got your random cell id.
x runs here from 0...nx and y from 0...ny.

Unfortunately this works only for regular shapes like yours, otherwise you get not only the error you see, but you also
would have to check for every cellID if the resulting space-point is inside the sensitive volume,
which makes it effectively unusable due to wasted CPU usage.

[ShujieL]

@MarkusFrankATcernch thanks! We did scan the cell ID and noticed that for a stave the segmentation occupies certain phase space, e.g. (0,6749) or (58786,65535) in seg_y (z direction). I am not sure though if it's a general solution b/c we do have Si disks formed with trapezoid, and moving forward will have inactive areas on each sensor unit, which may create discontinuity (not sure though, need your input).

Can you give me an example of how to check if a valid cell ID is indeed on a valid sensitive surface? We may not be able to afford it in production, but can at least try for initial test.

Eventually we need to generate electronic noise across all silicon detectors. are there alternative approaches we can explore?

[gaede]

The surfaces in DD4hep are attached to volumes that are used to define the surface boundaries. For drawing we the getLines() method that returns vectors of straight line on the surface boundary: https://github.com/AIDASoft/DD4hep/blob/master/DDRec/include/DDRec/Surface.h#L597-L600. This can obviously in cases of curved volumes not be exact, yet it should work for simple boxes or trapezoids. With these lines (or computing your own) it should be rather straight forward to sample uniform points on the surface inside the bounding lines (I believe) ...

[MarkusFrankATcernch]

@ShujieL I have no such example....
There is only a possibility in the Volumes to check if a point is inside. See TGeoVolume::Contains(...) ...

This is the method I think takes a lot if CPU.
If you have a trapezoid, you should be directly be able to compute the boundaries in terms of cell ids. This would be infinitely faster.