search

AIDASoft / DD4hep

Detector Description Toolkit for High Energy Physics
http://dd4hep.cern.ch
GNU Lesser General Public License v3.0
49 stars 95 forks source link

Tau decay products are produced at the origin with ddsim #1256

Closed Lrozzy closed 5 months ago

Lrozzy commented 5 months ago

Check duplicate issues.

Goal

Wanted to run an LLPs sensitivity study using susy-taus. Ran simulation over a healthy .hepmc file. The susy-taus decay to taus, but the tau decay products are being produced at the origin. Attached is a presentation showing more details/examples.

ddsim bug.pdf

Operating System and Version

Alma Linux 9

compiler

GCC 11

ROOT Version

6.28/02

DD4hep Version

1.25.1

Reproducer

I ran ddsim over the .hepmc file with the attached configurations and particle.tbl file.

The code I used to run ddsim is available at (specifically sim_Hbb/multi_sim.py): https://github.com/kdp-lab/Leo-LLPs-Code/tree/main

The particle.tbl, .hepmc, and steering (.py) files are available at: https://drive.google.com/drive/folders/1MjvZVrmsot-WBf74FLbVYZfWptAe-ZAX?usp=drive_link

Additional context

No response

andresailer commented 5 months ago

Is that from your pdf the case?

SIM.hepmc3.useHepMC3 = False

This doesn't work. Please use the hepmc3 reader.

Lrozzy commented 5 months ago

Is it ok to mix and match hepmc versions? For generation, the installation log in madgraph was hepmc2.06.09.

andresailer commented 5 months ago

Yes. HepMC3 has a hepmc2-ascii reader. The HepMC2 reader that is implement in DD4hep doesn't work for secondary vertices.

Lrozzy commented 5 months ago

Ok I see thank you. The reason we tried the hepmc2 reader is because using the hepmc3 reader, the taus were being produced at the origin. There is a screenshot in the presentation showing the taus coming from staus that decay far away from the origin but the taus are produced at the origin. We also do not see any stau or tau hits in simhit collections. Do you have any suggestions?

andresailer commented 5 months ago

This looks OK? The staus decay where the tau + Neutralino (?) start

[   id   ]index|      PDG |    px,     py,        pz    | px_ep,   py_ep , pz_ep      | energy  |gen|[simstat ]| vertex x,     y   ,   z     | endpoint x,    y  ,   z     |    mass |  charge |            spin             | colorflow | [parents] - [daughters]
[00053028]    4|  -2000015| 3.69e+03, 2.13e+03, 2.42e+03| 3.69e+03, 2.13e+03, 2.42e+03| 5.00e+03| 2 |[   t    ]| 0.00e+00, 0.00e+00, 0.00e+00| 1.16e+02, 6.73e+01, 7.63e+01| 1.00e+03| 1.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [1,3] - [6,7]
[00053029]    5|   2000015|-3.69e+03,-2.13e+03,-2.42e+03|-3.69e+03,-2.13e+03,-2.42e+03| 5.00e+03| 2 |[   t    ]| 0.00e+00, 0.00e+00, 0.00e+00|-1.66e+02,-9.63e+01,-1.09e+02| 1.00e+03|-1.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [1,3] - [8,9,24,25,26]
[00053030]    6|   1000049| 1.10e+03, 8.09e+02, 1.30e+03| 1.10e+03, 8.09e+02, 1.30e+03| 1.88e+03| 1 |[     l  ]| 1.16e+02, 6.73e+01, 7.63e+01| 2.54e+04, 1.87e+04, 3.00e+04| 1.00e-13| 0.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [4] - [] 
[00053031]    7|       -15| 2.59e+03, 1.32e+03, 1.12e+03| 2.59e+03, 1.32e+03, 1.12e+03| 3.12e+03| 2 |[    c   ]| 1.16e+02, 6.73e+01, 7.63e+01| 2.18e+02, 1.19e+02, 1.21e+02| 1.78e+00| 1.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [4] - [10,11,12,13]
[00053032]    8|   1000049|-2.78e+03,-1.61e+03,-2.29e+03|-2.78e+03,-1.61e+03,-2.29e+03| 3.95e+03| 1 |[     l  ]|-1.66e+02,-9.63e+01,-1.09e+02|-3.00e+04,-1.74e+04,-2.47e+04| 1.00e-13| 0.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [5] - []
[00053033]    9|        15|-9.07e+02,-5.18e+02,-1.30e+02|-9.07e+02,-5.18e+02,-1.30e+02| 1.05e+03| 2 |[    c   ]|-1.66e+02,-9.63e+01,-1.09e+02|-2.85e+02,-1.64e+02,-1.26e+02| 1.78e+00|-1.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [5] - [14,15,16]
[00053034]   10|       -16| 8.41e+02, 4.29e+02, 3.64e+02| 8.41e+02, 4.29e+02, 3.64e+02| 1.01e+03| 1 |[     l  ]| 2.18e+02, 1.19e+02, 1.21e+02| 3.00e+04, 1.53e+04, 1.30e+04| 0.00e+00| 0.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [7] - []
[00053035]   11|       111| 3.46e+02, 1.76e+02, 1.50e+02| 3.46e+02, 1.76e+02, 1.50e+02| 4.16e+02| 2 |[    c   ]| 2.18e+02, 1.19e+02, 1.21e+02| 2.18e+02, 1.19e+02, 1.21e+02| 1.35e-01| 0.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [7] - [17,18]
[00053036]   12|       111| 8.70e+02, 4.44e+02, 3.77e+02| 8.70e+02, 4.44e+02, 3.77e+02| 1.05e+03| 2 |[    c   ]| 2.18e+02, 1.19e+02, 1.21e+02| 2.18e+02, 1.19e+02, 1.21e+02| 1.35e-01| 0.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [7] - [19,20]
[00053037]   13|       211| 5.36e+02, 2.73e+02, 2.33e+02| 0.00e+00, 0.00e+00, 0.00e+00| 6.45e+02| 1 |[    c s ]| 2.18e+02, 1.19e+02, 1.21e+02| 1.13e+03, 5.82e+02, 5.15e+02| 1.40e-01| 1.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [7] - [40,41]
[00053038]   14|        16|-2.13e+02,-1.22e+02,-3.01e+01|-2.13e+02,-1.22e+02,-3.01e+01| 2.47e+02| 1 |[     l  ]|-2.85e+02,-1.64e+02,-1.26e+02|-3.00e+04,-1.72e+04,-4.32e+03| 0.00e+00| 0.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [9] - []
[00053039]   15|       311|-3.56e+02,-2.03e+02,-5.14e+01|-3.56e+02,-2.03e+02,-5.14e+01| 4.13e+02| 2 |[    c   ]|-2.85e+02,-1.64e+02,-1.26e+02|-2.85e+02,-1.64e+02,-1.26e+02| 4.98e-01| 0.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [9] - [21]

I used your 1000GeV_30mm.tbl file and set the GeneratorStatus of the staus from 22 to 2. ddsim doesn't pass anything with generator status > 2 to Geant4. What is important here also is the "SimStatus" which shows that the staus were actually simulated and decayed.

And I also get hits from the Staus, note the PDG of the MCParticle

--------------- print out of SimTrackerHit collection ---------------

  flag:  0x1d
 parameter CellIDEncoding [string]: system:0:8,barrel:8:3,layer:11:4,module:15:14,sensor:29:2,side:32:-2,strip:34:24,
     LCIO::THBIT_BARREL : 0
     LCIO::THBIT_MOMENTUM : 0
    quality bits: [os......]  o: hit from overlay s: hit from secondary not from the MCParticle associated to it

 [   id   ] |cellId0 |cellId1 |  position (x,y,z)               |   EDep   |   time   |PDG of MCParticle|   (px,  py, pz)   | pathLength  | Quality
------------|--------|--------|---------------------------------|----------|----------|-----------------|-------------------|-------------|---------
 [00054967] |537233409|00000000|(-2.37e+01, -1.37e+01, -1.55e+01)| 1.40e-05 | 1.07e-01 | 00000002000015  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,barrel:0,layer:1,module:11,sensor:1,side:0,strip:0)

 [00054968] |537235457|00000000|(-3.32e+01, -1.92e+01, -2.18e+01)| 1.93e-05 | 1.50e-01 | 00000002000015  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,barrel:0,layer:2,module:11,sensor:1,side:0,strip:0)

 [00054969] |537368577|00000000|(-4.43e+01, -2.56e+01, -2.91e+01)| 2.37e-05 | 2.01e-01 | 00000002000015  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,barrel:0,layer:3,module:15,sensor:1,side:0,strip:0)

 [00054970] |537468929|00000000|(-5.63e+01, -3.26e+01, -3.70e+01)| 2.05e-05 | 2.55e-01 | 00000002000015  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,barrel:0,layer:4,module:18,sensor:1,side:0,strip:0)

 [00054971] |537602049|00000000|(-6.75e+01, -3.90e+01, -4.43e+01)| 1.50e-05 | 3.05e-01 | 00000002000015  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,barrel:0,layer:5,module:22,sensor:1,side:0,strip:0)

 [00054972] |536938497|00000000|(+2.27e+01, +1.31e+01, +1.49e+01)| 1.82e-05 | 1.03e-01 | 000000-2000015  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,barrel:0,layer:1,module:2,sensor:1,side:0,strip:0)

 [00054973] |536940545|00000000|(+3.22e+01, +1.86e+01, +2.12e+01)| 1.54e-05 | 1.46e-01 | 000000-2000015  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,barrel:0,layer:2,module:2,sensor:1,side:0,strip:0)

 [00054974] |536975361|00000000|(+4.43e+01, +2.56e+01, +2.91e+01)| 1.11e-05 | 2.01e-01 | 000000-2000015  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,barrel:0,layer:3,module:3,sensor:1,side:0,strip:0)

 [00054975] |536977409|00000000|(+5.56e+01, +3.22e+01, +3.65e+01)| 1.23e-05 | 2.52e-01 | 000000-2000015  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,barrel:0,layer:4,module:3,sensor:1,side:0,strip:0)

 [00054976] |537012225|00000000|(+6.75e+01, +3.90e+01, +4.43e+01)| 1.31e-05 | 3.05e-01 | 000000-2000015  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,barrel:0,layer:5,module:4,sensor:1,side:0,strip:0)

I used ddsim from LCG_105b , so DD4hep version 1.27.02

Lrozzy commented 5 months ago

This looks great, thank you for trying it out. How did you change the GeneratorStatus of the staus?

andresailer commented 5 months ago

I used emacs to replace the 22 with 2 in the first event. Sed might also work

sed "s/\(2000015.* \)22 /\1 2 /" 1000_0.1.hepmc

This issue with the generator status showed up before, maybe we should add some treatment to allow this in ddsim.

Lrozzy commented 5 months ago

I just had a couple more questions:

  1. Can I still use dd4hep v1.25 for this or should I switch to 1.27?
  2. How did you get those printouts? Is that a script or something like an anajob that I can run in the command line?
andresailer commented 5 months ago
  1. I didn't see anything in the release notes that should prevent this from working with 1.25. When that issue with the status code has a solution it should be time to switch.
  2. These printouts comes from the lcio dumpevent command
andresailer commented 5 months ago

Hi @Lrozzy with the pull request #1260 you should be able to set

SIM.physics.alternativeDecayStatus = {22}

And run the unmodified hepmc file

andresailer commented 5 months ago

Also @Lrozzy , can you think of something to improve the documentation for this option (or other options) that would have made you realise this is what you need to set if it had already existed?

Lrozzy commented 5 months ago

Hi @Lrozzy with the pull request #1260 you should be able to set

SIM.physics.alternativeDecayStatus = {22}

And run the unmodified hepmc file

Thanks Andre. Unfortunately I think the status code 22 was only one of a few issues. Some of our staus promptly emit a photon and then continue on treated as different particles with different status codes:

Stau - unique id, status, endpoint: 743967 2 [0.0, 0.0]
Stau - unique id, status, endpoint: 743968 2 [0.0, 0.0]
Stau - unique id, status, endpoint: 743969 52 [180.24129240634133, -172.91154811165617]
Stau - unique id, status, endpoint: 743970 51 [11.41540434665874, 10.95101481738347]

In this case the status of the staus that continue on are 51 and 52. These are defined as:

When we just changed the status code 22 --> 2, we realised that we also have taus that come from staus but, again, promptly emit a photon and continue on, treated as different particles with different status codes:

Tau Parent - unique id, PDG, endpoint: [(743969, -2000015, [-99.03976248015312, -150.59232694989174, -172.91154811165617])] 
Tau - unique id, production vertex, endpoint, status, : 743973 [180.24129240634133, -172.91154811165617] [180.24129240634133, -172.91154811165617] 23 
        Tau decay products - unique id, PDG, status, prod_vertex:
         (743977, -15, 2, [0.0, 0.0, 0.0])
         (743978, 22, 1, [0.0, 0.0, 0.0])

Tau Parent - unique id, PDG, endpoint: [(743973, -15, [180.24129240634133, -172.91154811165617])] 
Tau - unique id, vertex, endpoint, status, : 743977 [0.0, 0.0] [38.13097318884068, -38.02263542816206] 2
                 Tau decay products - unique id, PDG, status, prod_vertex:
         (743982, -16, 1, [-14.262259419882717, -35.363244655546715, -38.02263542816206])
         (743983, 111, 2, [-14.262259419882717, -35.363244655546715, -38.02263542816206])
         (743984, 211, 1, [-14.262259419882717, -35.363244655546715, -38.02263542816206])
         (743985, 311, 2, [-14.262259419882717, -35.363244655546715, -38.02263542816206])

In this case you can see the tau that emits a photon has status code 23, and the one that continues has status code 2, but we do not see the tau in simhits. Also, somehow the tau that promptly emits a photon emits it at the origin despite being very far away from the origin, but this is not a problem for the status 2 tau that continues on.

andresailer commented 5 months ago

HI @Lrozzy Can you give me the event number for this? SIM.physics.alternativeDecayStatus = {22, 51, 52,...} might work then.

Lrozzy commented 5 months ago

Sure, this was event number 1, index 1. So the 2nd event of the 1000_0.1.hepmc file

andresailer commented 5 months ago

@Lrozzy 

SIM.physics.alternativeDecayStatuses = {22, 23, 51, 52}
SIM.physics.zeroTimePDGs = {17, 11} # removing 13, 15 from this list

Seems to work for me with event 1

Lrozzy commented 5 months ago

Ok great, I'll try and do the same!

andresailer commented 5 months ago

If the applied fixes are not enough, please let us know. (Also let us know if everything works).

Lrozzy commented 4 months ago

Just a quick update - we got ddsim working and the output seems ok (stau and tau hits present) but when I try to run digi it fails, using both the old stack and the new. I've attached the outputs for both. Please let me know if you have any suggestions or want more information.

old_stack_output.log new_output.log

andresailer commented 4 months ago

Hi @Lrozzy Great that the simulation works now!

The digi issue however is outside the scope for the DD4hep repository.

Just some comments: Old:

      /tmp/root/spack-stage/spack-stage-lcio-2.19.1-3sauiiqpg4ch56alksvwg46wvqb6aele/spack-src/src/cpp/src/SIO/SIOWriter.cc (l.91) in open: [SIOWriter::open()] Couldn't open file: '/local/d1/mu+mu-/digi/4000_0.1_digi.slcio' [not_open]

File not found? Can't help you there.

New:

Warning in <TClassTable::Add>: class _Rb_tree_iterator<pair<const string,map<string,string> > > already in TClassTable
...

Points towards a mixture of environments: same packages with different versions... So you need to rebuild your whole stack with the new DD4hep version in a consistent way?

As I don't know who or how you build the stack, I can't tell you where would be a better place to raise this issue.

Cheers, Andre

Lrozzy commented 4 months ago

Thanks @andresailer. The file not found is strange because when I try to run digi on a sim file produced from the old stack (dd4hep 1.25) it works fine, and there shouldn't be any permission issues. As for the new, I built dd4hep and lcgeo, and then I sourced first the cvmfs/muoncollider, then my dd4hep and my lcgeo. This was how I made the working sim file.

source /cvmfs/muoncollider.cern.ch/release/2.8-patch2/setup.sh
source /local/d1/lrozanov/mucoll-tutorial-2023/DD4hep/bin/thisdd4hep.sh
source /local/d1/lrozanov/mucoll-tutorial-2023/lcgeo/bin/thislcgeo.sh
andresailer commented 4 months ago 
source /cvmfs/muoncollider.cern.ch/release/2.8-patch2/setup.sh
source /local/d1/lrozanov/mucoll-tutorial-2023/DD4hep/bin/thisdd4hep.sh
source /local/d1/lrozanov/mucoll-tutorial-2023/lcgeo/bin/thislcgeo.sh

This is exactly what I mean with mixing environments. It works for simulation, because you are just using DD4hep, but your Digi was build against an older DD4hep version, so you start mixing libraries.

Please open a new issue with all the details that would allow someone to reproduce what your are doing.