Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
1 task done
Lrozzy opened this issue Apr 29, 2024 · 21 comments · Fixed by #1260
Closed
1 task done

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

Lrozzy opened this issue Apr 29, 2024 · 21 comments · Fixed by #1260
Assignees

Comments

@Lrozzy
Copy link

Lrozzy commented Apr 29, 2024

Check duplicate issues.

  • Checked for duplicates

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

@Lrozzy Lrozzy added the bug label Apr 29, 2024
@andresailer
Copy link
Member

Is that from your pdf the case?

SIM.hepmc3.useHepMC3 = False

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

@Lrozzy
Copy link
Author

Lrozzy commented Apr 29, 2024

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

@andresailer
Copy link
Member

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

@Lrozzy
Copy link
Author

Lrozzy commented Apr 29, 2024

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
Copy link
Member

  • I am trying to get the staus to "fly"...

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
Copy link
Author

Lrozzy commented Apr 29, 2024

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

@andresailer
Copy link
Member

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
Copy link
Author

Lrozzy commented Apr 30, 2024

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
Copy link
Member

  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
Copy link
Member

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

SIM.physics.alternativeDecayStatus = {22}

And run the unmodified hepmc file

@andresailer
Copy link
Member

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
Copy link
Author

Lrozzy commented May 7, 2024

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:

  • 51 : outgoing produced by parton branching
  • 52 : outgoing copy of recoiler, with changed momentum

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
Copy link
Member

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

@Lrozzy
Copy link
Author

Lrozzy commented May 7, 2024

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

@andresailer
Copy link
Member

@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
Copy link
Author

Lrozzy commented May 8, 2024

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

@andresailer
Copy link
Member

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

@Lrozzy
Copy link
Author

Lrozzy commented May 24, 2024

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
Copy link
Member

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
Copy link
Author

Lrozzy commented May 28, 2024

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
Copy link
Member

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants