Merge pull request cms-sw#29 from ajgilbert/add-pythia8-gun

Add option to use pythia8 particle gun

README.md

@@ -26,10 +26,11 @@ To produce `NEVENTS` GEN-SIM-DIGI events with `NPART` sets of particles (per eve
   --thresholdMin PTMIN
   --thresholdMax PTMAX
   --gunType GUNTYPE
+  [--gunMode GUNMODE]
   [--local]
   --tag MYTAG
 ```
-Here, one can produce a custom set of particles by providing `PART_PDGID` as a set of comma-separated single PDG IDs.
+Here, one can produce a custom set of particles by providing `PART_PDGID` as a set of comma-separated single PDG IDs. To simulate the decay of unstable particles, e.g. quarks, gluons or taus,  an alternative particle gun based on PYTHIA8 can be used by setting `--gunMode pythia8`. 
 
 To produce `NEVENTS` GEN-SIM-DIGI events with pair of particles within given angular distance ΔR(η,φ) (per event), where the first particle is of type `PART_PDGID` and in the p_T range from `PTMIN` to `PTMAX`, and the second one is of type `INCONE_PART_PDGID` and at distance from `DRMIN` to `DRMAX` and with p_T in range from `PTRATIO_MIN` to `PTRATIO_MAX` relative to the first particle, one should run:
 ```
@@ -43,6 +44,7 @@ To produce `NEVENTS` GEN-SIM-DIGI events with pair of particles within given ang
   --thresholdMin PTMIN
   --thresholdMax PTMAX
   --gunType Pt
+  [--gunMode GUNMODE]
   --InConeID INCONE_PART_PDGID
   --MinDeltaR DRMIN
   --MaxDeltaR DRMAX
@@ -51,7 +53,7 @@ To produce `NEVENTS` GEN-SIM-DIGI events with pair of particles within given ang
   [--local]
   --tag MYTAG
 ```
-One should also note that for the genertion of pairs of particles within a given cone, one has to use the "Pt" gun.
+One should also note that for the genertion of pairs of particles within a given cone, one has to use the "Pt" gun. Also note that it is currently not possible to generate pairs within a cone using the PYTHIA8-based gun.
 
 The script will create a directory called `partGun_[MYTAG]_[DATE]` locally or on the CMG EOS area (see options), and submit jobs to queue `QUEUENAME` with `NPERJOB` events per job,
 `NEVENTS` in total.

SubmitHGCalPGun.py

@@ -28,6 +28,7 @@ def parseOptions():
     parser.add_option('', '--nPart',  dest='NPART',  type=int,   default=1,      help='number of times particles of type(s) PARTID will be generated per event, default is 1')
     parser.add_option('', '--thresholdMin',  dest='thresholdMin',  type=float, default=1.0,     help='min. threshold value')
     parser.add_option('', '--thresholdMax',  dest='thresholdMax',  type=float, default=35.0,    help='max. threshold value')
+    parser.add_option('', '--gunMode',   dest='gunMode',   type='string', default='default',    help='default or pythia8')
     parser.add_option('', '--gunType',   dest='gunType',   type='string', default='Pt',    help='Pt or E gun')
     parser.add_option('', '--PU',   dest='PU',   type='string', default='0',    help='PU value (0 is the default)')
     parser.add_option('', '--PUDS',   dest='PUDS',   type='string', default='',    help='PU dataset')
@@ -59,6 +60,11 @@ def parseOptions():
         print 'ERROR: CMSSW does not seem to be set up. Exiting...'
         sys.exit()
 
+    partGunModes = ['default', 'pythia8']
+    if opt.gunMode not in partGunModes:
+        parser.error('Particle gun mode ' + opt.gunMode + ' is not supported. Exiting...')
+        sys.exit()
+
     partGunTypes = ['Pt', 'E']
     if opt.gunType not in partGunTypes:
         parser.error('Particle gun type ' + opt.gunType + ' is not supported. Exiting...')
@@ -80,7 +86,7 @@ def parseOptions():
 
     # list of supported particles, check if requested partID list is a subset of the list of the supported ones
     global particles
-    particles = ['22', '111', '211', '11', '13', '15', '12', '14', '16', '130']
+    particles = ['22', '111', '211', '11', '13', '15', '12', '14', '16', '130', '1', '2', '3', '4', '5', '22']
     inPartID = [p.strip(" ") for p in opt.PARTID.split(",")] # prepare list of requested IDs (split by ",", strip white spaces)
     if not (set(inPartID) < set(particles) or opt.PARTID == ''):
         parser.error('Particle(s) with ID(s) ' + opt.PARTID + ' is not supported. Exiting...')
@@ -137,7 +143,7 @@ def printSetup(CMSSW_BASE, CMSSW_VERSION, SCRAM_ARCH, currentDir, outDir):
         curr_input= opt.RELVAL
     print 'PU:         ',opt.PU
     print 'PU dataset: ',opt.PUDS
-    print 'INPUTS:     ', [curr_input, 'Particle gun type: ' + opt.gunType + ', PDG ID '+str(opt.PARTID)+', '+str(opt.NPART)+' times per event, ' + opt.gunType + ' threshold in ['+str(opt.thresholdMin)+','+str(opt.thresholdMax)+']',opt.RELVAL][int(opt.DTIER=='GSD')]
+    print 'INPUTS:     ', [curr_input, 'Particle gun mode: ' + opt.gunMode + ', type: ' + opt.gunType + ', PDG ID '+str(opt.PARTID)+', '+str(opt.NPART)+' times per event, ' + opt.gunType + ' threshold in ['+str(opt.thresholdMin)+','+str(opt.thresholdMax)+']',opt.RELVAL][int(opt.DTIER=='GSD')]
     if (opt.InConeID!='' and opt.DTIER=='GSD'):
         print '             IN-CONE: PDG ID '+str(opt.InConeID)+', deltaR in ['+str(opt.MinDeltaR)+ ','+str(opt.MaxDeltaR)+']'+', momentum ratio in ['+str(opt.MinMomRatio)+ ','+str(opt.MaxMomRatio)+']'
     print 'STORE AREA: ', [opt.eosArea, currentDir][int(opt.LOCAL)]
@@ -183,6 +189,8 @@ def submitHGCalProduction():
     SCRAM_ARCH = os.getenv('SCRAM_ARCH')
     commonFileNamePrefix = 'partGun'
     partGunType = 'FlatRandom%sGunProducer' % opt.gunType
+    if opt.gunMode == 'pythia8':
+        partGunType = 'Pythia8%sGun' % opt.gunType
     if opt.InConeID != '':
         partGunType = 'MultiParticleInConeGunProducer'  # change part gun type if needed, keep opt.gunType unchanged (E or Pt) for the "primary particle"
 
@@ -361,6 +369,7 @@ def submitHGCalProduction():
                 s_template=s_template.replace('MAXTHRESHSTRING',"Max"+str(opt.gunType))
                 s_template=s_template.replace('MINTHRESHSTRING',"Min"+str(opt.gunType))
                 s_template=s_template.replace('DUMMYPU',str(mixing))
+                s_template=s_template.replace('GUNMODE',str(opt.gunMode))
 
 
             elif (opt.DTIER == 'RECO' or opt.DTIER == 'NTUP'):

templates/partGun_GSD_template.py

@@ -78,22 +78,42 @@
 from Configuration.AlCa.GlobalTag import GlobalTag
 process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic', '')
 
-process.generator = cms.EDProducer("GUNPRODUCERTYPE",
-    AddAntiParticle = cms.bool(True),
-    PGunParameters = cms.PSet(
-        MaxEta = cms.double(3.0),
-        MaxPhi = cms.double(3.14159265359),
-        MAXTHRESHSTRING = cms.double(DUMMYTHRESHMAX),
-        MinEta = cms.double(1.479),
-        MinPhi = cms.double(-3.14159265359),
-        MINTHRESHSTRING = cms.double(DUMMYTHRESHMIN),
-        #DUMMYINCONESECTION
-        PartID = cms.vint32(DUMMYIDs)
-    ),
-    Verbosity = cms.untracked.int32(0),
-    firstRun = cms.untracked.uint32(1),
-    psethack = cms.string('multiple particles predefined pT/E eta 1p479 to 3')
-)
+gunmode = 'GUNMODE'
+
+if gunmode == 'default':
+    process.generator = cms.EDProducer("GUNPRODUCERTYPE",
+        AddAntiParticle = cms.bool(True),
+        PGunParameters = cms.PSet(
+            MaxEta = cms.double(3.0),
+            MaxPhi = cms.double(3.14159265359),
+            MAXTHRESHSTRING = cms.double(DUMMYTHRESHMAX),
+            MinEta = cms.double(1.479),
+            MinPhi = cms.double(-3.14159265359),
+            MINTHRESHSTRING = cms.double(DUMMYTHRESHMIN),
+            #DUMMYINCONESECTION
+            PartID = cms.vint32(DUMMYIDs)
+        ),
+        Verbosity = cms.untracked.int32(0),
+        firstRun = cms.untracked.uint32(1),
+        psethack = cms.string('multiple particles predefined pT/E eta 1p479 to 3')
+    )
+elif gunmode == 'pythia8':
+    process.generator = cms.EDFilter("GUNPRODUCERTYPE",
+        maxEventsToPrint = cms.untracked.int32(1),
+        pythiaPylistVerbosity = cms.untracked.int32(1),
+        pythiaHepMCVerbosity = cms.untracked.bool(True),
+        PGunParameters = cms.PSet(
+          ParticleID = cms.vint32(DUMMYIDs),
+          AddAntiParticle = cms.bool(True),
+          MinPhi = cms.double(-3.14159265359),
+          MaxPhi = cms.double(3.14159265359),
+          MINTHRESHSTRING = cms.double(DUMMYTHRESHMIN),
+          MAXTHRESHSTRING = cms.double(DUMMYTHRESHMAX),
+          MinEta = cms.double(1.479),
+          MaxEta = cms.double(3.0)
+          ),
+        PythiaParameters = cms.PSet(parameterSets = cms.vstring())
+    )
 
 # Path and EndPath definitions
 process.generation_step = cms.Path(process.pgen)