DCH v2 digitizer, which smears the position and adds cluster information

[atolosadelgado]

BEGINRELEASENOTES

• New digitizer for DCH v2 (based on twisted tubes). It calculates cluster information and smears the position along and perpendicular the wire separately. New custom EDM4hep data extension was added to include these data.

ENDRELEASENOTES

Hi,

Please find attached the drift chamber (DCH) digitizer, based on Walaa's code for the cluster number and size calculations. In addition, I added the position smearing based on analytical calculations (I was stuck here, but thanks to Andre we found the problem).

The algorithm works as expected, but I can imagine there is room for improvement regarding naming, coding style in general, so please let me know if I can improve it somehow.

I have been working on a stand alone repository, in case you want to try out: https://github.com/atolosadelgado/DCH_detector

Alvaro.

[BrieucF]

Hi Alvaro, excellent! Can you add a test and add a link between sim and digi hits? See e.g. https://github.com/key4hep/k4RecTracker/pull/28/files

[atolosadelgado]

Hi Alvaro, excellent! Can you add a test and add a link between sim and digi hits? See e.g. https://github.com/key4hep/k4RecTracker/pull/28/files

I have added the association and the test. I have inspected the output file and the associations seems to make sense.

Associations are now an additional output collection of the algorithm, but I am not sure if this is the proper way.

You can use the following script to inspect the output file of the test, everything seems consistent

import ROOT
from podio import root_io
import dd4hep as dd4hepModule
from ROOT import dd4hep
from ROOT import TVector3, TFile, TH1D

ofile = TFile("kk.root","update")
hDpw_calc = TH1D("hDpw_calc","",100,0,10)
hDpw_alg  = TH1D("hDpw_alg" ,"",100,0,10)
hDpw_diff = TH1D("hDpw_diff","",100,0,10)

input_file_path = "dch_proton_10GeV_digi.root"
podio_reader = root_io.Reader(input_file_path)

for event in podio_reader.get("events"):
  for dc_hit in event.get("DCH_DigiSimAssociationCollection"):
    digihit =  dc_hit.getDigi()
    simhit  =  dc_hit.getSim()
    digipos = TVector3( digihit.getPosition()[0],digihit.getPosition()[1],digihit.getPosition()[2] )
    simpos  = TVector3(  simhit.getPosition()[0], simhit.getPosition()[1], simhit.getPosition()[2] )
    hit_wire_vector = simpos - digipos
    distance_calculated = hit_wire_vector.Mag()
    distance_from_alg = digihit.getDistanceToWire()
    # print( distance_from_alg - distance_calculated )
    #print( f"{digihit.getEDep()}      {simhit.getEDep()}")
    hDpw_calc.Fill(distance_calculated)
    hDpw_alg.Fill(distance_from_alg)
    hDpw_diff.Fill(distance_calculated-distance_from_alg)
    
ofile.cd()
hDpw_calc.Write()
hDpw_alg.Write()
hDpw_diff.Write()
ofile.Close()

[atolosadelgado]

@BrieucF the code is ready for review :)

[BrieucF]

Thanks Alvaro!

[BrieucF]

Let's add the interface as well: https://github.com/BrieucF/k4RecTracker/blob/fix_for_Dolores/DCHdigi/dataFormatExtension/driftChamberHit.yaml#L66

[atolosadelgado]

Let's add the interface as well: https://github.com/BrieucF/k4RecTracker/blob/fix_for_Dolores/DCHdigi/dataFormatExtension/driftChamberHit.yaml#L66

does this do the trick? how can I test that it works?
47b0d5d

[BrieucF]

Let's add the interface as well: https://github.com/BrieucF/k4RecTracker/blob/fix_for_Dolores/DCHdigi/dataFormatExtension/driftChamberHit.yaml#L66

does this do the trick? how can I test that it works? 47b0d5d

I am surprised that it does not complain that the interfaces edm4hep::TrackerHit: is then defined twice in edm4hep. Let's ask @tmadlener what he thinks. Context: the goal is to avoid having to redefine extension::TrackerHit3D and extension::TrackerHitPlane in order to build an interface, as done here https://github.com/BrieucF/k4RecTracker/blob/fix_for_Dolores/DCHdigi/dataFormatExtension/driftChamberHit.yaml#L66. Last time I tried I could not mix extension objects and edm4hep object in an interface. The latter is indeed very annoying because, in order for subsequent algorithm to work, we have to cast e.g. edm4hep::TrackerHit3D coming from the SIM step into extension::TrackerHit3D and store this duplicated collection in the output rootfile.

[BrieucF]

All these parameters will be set at every call of the function, since many are hard coded, would'nt it be better to define them in the header once and for all?

[atolosadelgado]

Let's add the interface as well: https://github.com/BrieucF/k4RecTracker/blob/fix_for_Dolores/DCHdigi/dataFormatExtension/driftChamberHit.yaml#L66

does this do the trick? how can I test that it works? 47b0d5d

I am surprised that it does not complain that the interfaces edm4hep::TrackerHit: is then defined twice in edm4hep. Let's ask @tmadlener what he thinks. Context: the goal is to avoid having to redefine extension::TrackerHit3D and extension::TrackerHitPlane in order to build an interface, as done here https://github.com/BrieucF/k4RecTracker/blob/fix_for_Dolores/DCHdigi/dataFormatExtension/driftChamberHit.yaml#L66. Last time I tried I could not mix extension objects and edm4hep object in an interface. The latter is indeed very annoying because, in order for subsequent algorithm to work, we have to cast e.g. edm4hep::TrackerHit3D coming from the SIM step into extension::TrackerHit3D and store this duplicated collection in the output rootfile.

I changed the data extension as you proposed originally (by duplicating code), so the interface can work as expected.

[atolosadelgado]

Hi,

I have changed the data extension for this digitizer, so the size of each cluster is given (before the total number of electrons along the step was reported)

Few remarks:

• It makes a smearing in the position along the wire and radially, according to the input parameters
• Number of clusters and their size is calculated as the option 3 from this paper [1]
• Debug histograms are created if debug option is enabled
• Stand alone test run simulation of the drift chamber based on twisted tubes, and then apply the digitizer. Dedicated directory with all the files needed is given in DCHdigi/test/test_DCHdigi/
• It requires that the cellID contain the layer and number of cell within the layer (nphi). It does not matter if the segmentation comes from geometrical segmentation by using twisted tubes and hyperboloids (and the cellID is created out of volume IDs), or the segmentation is virtual DD4hep segmentation.

Some plots:

1. individual cluster size for kaons at 45 GeV
   
2. individual cluster size for protons at 5 GeV
   
3. impact of step length, comparison 1mm vs 10 mm (default range cut 0.7mm, for both), protons at 10 GeV with direction (1,1,1), discrepancies on 1-5 electron-size cluster number, otherwise shape is similar.
   
4. distance of hit to the wire (protons 10GeV, direction 111, max step length 1 mm)
   
5. Comparison of Z-position of Digi (projected + smeared) hit vs Sim hit. Main contribution to the change is that the position of the digitized hit is the projection of the sim hit position onto the wire.

TFile *_file0 = TFile::Open("dch_proton_10GeV_digi.root")
events->Draw("DCH_DigiCollection.position.z - DCHCollection.position.z")

5.1 Comparison of Z position of digitized hits vs simhits, no smearing applied (just projection)

5.2 Comparison, enabling smearing along the wire (1mm) and perpendicularly (0.1mm)

[atolosadelgado]

I have added a suffix with the version v01 to this new digitizer, and writen a readme file in the DCHdigi directory

please let me know if there is anything else I can do :)

[BrieucF]

Thanks a lot Alvaro! Can you add a plot of "distance to the wire of digitized hit - distance of the wire from the simhit" together with the code snippet that allows to derive the distance to the wire from the simHit position + cellID (using the facilities from DCH_info)? By the way, can this also be done in Python?

[atolosadelgado]

distance to the wire of digitized hit

Hi,

Thanks for reviewing it!

I am not sure I understood your question:

• distance to the wire of digitized hit is zero, by construction: the digitized hit position is a point along the closest wire to the sim hit.
• distance of the wire from the simhit, this information is stored in the output class in the member named distanceToWire, line 58 of the yaml file
• To draw a 2D histogram with the distance to the wire and the cell ID, we can use:

events->Draw("DCH_DigiCollection.distanceToWire:DCH_DigiCollection.cellID")

For more advanced plots, I guess one would need a short script, which can be in python

[BrieucF]

distance to the wire of digitized hit

Hi,

Thanks for reviewing it!

I am not sure I understood your question:

* `distance to the wire of digitized hit` is zero, by construction: the digitized hit position is a point along the closest wire to the sim hit.

What is "zero" is the distance between the position member of the digitizedHit and the wire, by definition. What I am referring to here is given by the member distanceToWire of the digitized hit

* `distance of the wire from the simhit`, this information is stored in the output class in the member named `distanceToWire`, [line 58 of the yaml file](https://github.com/key4hep/k4RecTracker/blob/39555a8792f73ba5ad93b0b736a061c72e0cc840/DCHdigi/dataFormatExtension/driftChamberHit.yaml#L58)

No, the information you point is for the digitized hit. For the simHit, we need a routine (I thought you had one defined in the DCH_info) that gives the distance to the wire using the CellID and the simHit position.

* To draw a 2D histogram with the distance to the wire and the cell ID, we can use:

events->Draw("DCH_DigiCollection.distanceToWire:DCH_DigiCollection.cellID")

For more advanced plots, I guess one would need a short script, which can be in python

The plot I am asking for is a validation: getDistanceToTheWire(simHit.position, cellID) - digitizedHit.getDistanceToWire. And we should see a gaussian with 100 um. Indirectly this code snippet will also provide an example on how to retrieve this info for simHit.

Is it clearer?

[BrieucF]

Attaching here the number of cluster for 10 GeV protons as produced by the tests

[BrieucF]

Thanks Alvaro!