JeffersonLab / hps-mc

HPS MC toolkit
1 stars 6 forks source link

add iDM particles to displacement program #387

Closed tomeichlersmith closed 1 year ago

tomeichlersmith commented 1 year ago

The iDM model I've developed has different PDG numbers compared to the previous dark photon and SIMP models so a small edit needs to be made to the displacement program. Besides this inclusion, I also did a very small amount of code cleanup: removing unused variables and updating comments.

Evidence

I ran hps-mc-job run -d working idm_job.py job.json with this branch using

idm_job.py

from hpsmc.generators import MG5
from hpsmc.tools import SLIC, Unzip, DisplaceUni, BeamCoords, FilterBunches, JobManager

job.description = 'iDM generation'
job.add([
    MG5(name='idm', event_types=['unweighted']),
    Unzip(inputs=['idm_unweighted_events.lhe.gz'], outputs=['idm.lhe']),
    DisplaceUni(inputs=['idm.lhe'], outputs=['idm.stdhep'])
])

job.json

{
    "nevents": 10000,
    "run_params": "2pt3",
    "Map": 300.0,
    "mchi": 100.0,
    "dmchi": 10.0,
    "seed": 879483
}

and then I could print out the stdhep file and see the produced electrons have a displacement.

$ stdhep_print_stdhep working/idm.stdhep 1
Reading from idm.stdhep; expecting 1000000 events
 StdHepXdrReadOpen: successfully opened input stream 1
          title: idm.stdhep
          date: Thu May 25 09:51:32 2023

                    10002 events
                    8 blocks per event
 mcfio_Block: Unable to find block i.d. 101 in Stream 1 
 mcfio_Block: Unable to find block i.d. 105 in Stream 1 
ilbl = 100
read event 1: nevhep = 1, nhep = 10
isthep = 3      idhep = 11      jmohep[0] = 0   jmohep[1] = 0   jdahep[0] = 3   jdahep[1] = 6
phep[0] = 0.000000      phep[1] = 0.000000      phep[2] = 2.300000      phep[3] = 2.300000      phep[4] = 0.000511
vhep[0] = 0.000000      vhep[1] = 0.000000      vhep[2] = 0.000000      vhep[3] = 0.000000
isthep = 3      idhep = 9000002 jmohep[0] = 0   jmohep[1] = 0   jdahep[0] = 3   jdahep[1] = 6
phep[0] = -0.000000     phep[1] = -0.000000     phep[2] = -0.000000     phep[3] = 174.000000    phep[4] = 174.000000
vhep[0] = 0.000000      vhep[1] = 0.000000      vhep[2] = 0.000000      vhep[3] = 0.000000
isthep = 2      idhep = 1023    jmohep[0] = 1   jmohep[1] = 2   jdahep[0] = 4   jdahep[1] = 10
phep[0] = 0.021856      phep[1] = -0.015166     phep[2] = 2.225261      phep[3] = 2.245516      phep[4] = 0.299746
vhep[0] = 0.000000      vhep[1] = 0.000000      vhep[2] = 0.000000      vhep[3] = 0.000000
isthep = 2      idhep = 1000023 jmohep[0] = 3   jmohep[1] = 3   jdahep[0] = 7   jdahep[1] = 9
phep[0] = 0.082454      phep[1] = 0.051219      phep[2] = 1.579686      phep[3] = 1.586827      phep[4] = 0.114847
vhep[0] = 0.000000      vhep[1] = 0.000000      vhep[2] = 0.000000      vhep[3] = 0.000000
isthep = 1      idhep = 11      jmohep[0] = 1   jmohep[1] = 2   jdahep[0] = 0   jdahep[1] = 0
phep[0] = -0.031058     phep[1] = -0.001391     phep[2] = 0.044736      phep[3] = 0.054480      phep[4] = 0.000511
vhep[0] = 0.000000      vhep[1] = 0.000000      vhep[2] = 0.000000      vhep[3] = 0.000000
isthep = 1      idhep = 9000002 jmohep[0] = 1   jmohep[1] = 2   jdahep[0] = 0   jdahep[1] = 0
phep[0] = 0.009202      phep[1] = 0.016557      phep[2] = 0.030003      phep[3] = 174.000004    phep[4] = 174.000000
vhep[0] = 0.000000      vhep[1] = 0.000000      vhep[2] = 0.000000      vhep[3] = 0.000000
isthep = 1      idhep = 1000022 jmohep[0] = 4   jmohep[1] = 4   jdahep[0] = 0   jdahep[1] = 0
phep[0] = 0.082050      phep[1] = 0.045592      phep[2] = 1.312819      phep[3] = 1.319594      phep[4] = 0.095000
vhep[0] = 3.258894      vhep[1] = 2.024387      vhep[2] = 62.435541     vhep[3] = 62.717779
isthep = 1      idhep = -11     jmohep[0] = 4   jmohep[1] = 4   jdahep[0] = 0   jdahep[1] = 0
phep[0] = 0.000499      phep[1] = -0.003980     phep[2] = 0.106805      phep[3] = 0.106881      phep[4] = 0.000511
vhep[0] = 3.258894      vhep[1] = 2.024387      vhep[2] = 62.435541     vhep[3] = 62.717779
isthep = 1      idhep = 11      jmohep[0] = 4   jmohep[1] = 4   jdahep[0] = 0   jdahep[1] = 0
phep[0] = -0.000096     phep[1] = 0.009607      phep[2] = 0.160062      phep[3] = 0.160351      phep[4] = 0.000511
vhep[0] = 3.258894      vhep[1] = 2.024387      vhep[2] = 62.435541     vhep[3] = 62.717779
isthep = 1      idhep = 1000022 jmohep[0] = 3   jmohep[1] = 3   jdahep[0] = 0   jdahep[1] = 0
phep[0] = -0.060598     phep[1] = -0.066386     phep[2] = 0.645575      phep[3] = 0.658689      phep[4] = 0.095000
vhep[0] = 0.000000      vhep[1] = 0.000000      vhep[2] = 0.000000      vhep[3] = 0.000000

Lack of Side Effects

This will not affect nominal A' and SIMP production since neither of those models have a particle with the PDG that iDM requires so simply including it in an OR clause will support both production methods without effecting those.

Now, there are models with a PDG matching this particle - the MSSM models have a neutralino with this PDG - so if supporting displacement of these models with this program is necessary I would need to update it.

tomeichlersmith commented 1 year ago

Yes, all of the new files are a result of PR #366 on which this is based, so any discussion about how we are storing the MG generation type should be done there.

cbravo135 commented 1 year ago

I think this is a good PR to request you to add idm_job.py to python/jobs and make an idm example in the examples directory

tomeichlersmith commented 1 year ago

I have been able to confirm that this displacement is picked up by slic. The following plot is draw from the histogram filled by hpstr in the hps-mc job below. Since both electron and positron have the same vertex, only one is visible.

vertex_z

idm_job.py

from hpsmc.generators import MG5
from hpsmc.tools import SLIC, Unzip, DisplaceUni, BeamCoords, FilterBunches, JobManager

job.description = 'iDM generation'
job.add([
    MG5(name='idm', event_types=['unweighted']),
    Unzip(inputs=['idm_unweighted_events.lhe.gz'], outputs=['idm.lhe']),
    DisplaceUni(inputs=['idm.lhe'], outputs=['idm.stdhep']),
    BeamCoords(),
    SLIC(),
    HPSTR(cfg = 'tuplize', inputs='idm_rot.slcio', outputs='idm_rot.root'),
    HPSTR(cfg = 'fill', inputs='idm_rot.root', outputs='h_idm_rot.root'),
    FilterBunches(),
    JobManager(steering='readout'),
    JobManager(steering='recon')
])

job.json

{
    "nevents": 10000,
    "run_params": "2pt3",
    "Map": 300.0,
    "mchi": 100.0,
    "dmchi": 10.0,
    "seed": 879483,
    "target_thickness": 0.002,
    "beam_energy": 2300.0,
    "num_electrons": 1500,
    "target_z": -4.3,
    "run_number": 7800,
    "detector": "HPS-PhysicsRun2016-Pass2",
    "steering_files": {
      "readout": "/org/hps/steering/readout/PhysicsRun2016TrigPairs1.lcsim",
      "recon": "/org/hps/steering/recon/PhysicsRun2016FullReconMC.lcsim"
    },
    "config_files": {
      "tuplize" : "mcTuple_cfg.py",
      "fill"    : "anaMCTuple_cfg.py"
    }
}