File indexing completed on 2024-10-03 05:27:22
0001
0002 '''
0003 This script is designed to read the G4msg.log output of a cmsDriver.py command run in a test release with a special configuration of SimG4Core/Application/RunManager.cc and a special version of SimulationG4.py fragment.
0004 To dump the G4 particle table one needs to:
0005 1-Edit SimG4Core/Application/src/RunManager.cc:
0006 a-Add in the includes at the top:
0007 #include "G4ParticleTable.hh"
0008 b-Add at the bottom of the RunManager::initG4 method:
0009 edm::LogInfo("SimG4CoreApplication") << "Output of G4ParticleTable DumpTable:";
0010 G4ParticleTable::GetParticleTable()->DumpTable("ALL");
0011 2-Edit the Validation/Performance/python/TimeMemoryG4Info.py customise fragment (or you could create your own):
0012 a-Configure the output (in this case to the file G4msg.log) to include SimG4CoreApplication:
0013 process.MessageLogger.files = dict(G4msg = cms.untracked.PSet(
0014 noTimeStamps = cms.untracked.bool(True)
0015 #First eliminate unneeded output
0016 ,threshold = cms.untracked.string('INFO')
0017 ,INFO = cms.untracked.PSet(limit = cms.untracked.int32(0))
0018 ,FwkReport = cms.untracked.PSet(limit = cms.untracked.int32(0))
0019 ,FwkSummary = cms.untracked.PSet(limit = cms.untracked.int32(0))
0020 ,Root_NoDictionary = cms.untracked.PSet(limit = cms.untracked.int32(0))
0021 ,FwkJob = cms.untracked.PSet(limit = cms.untracked.int32(0))
0022 ,TimeReport = cms.untracked.PSet(limit = cms.untracked.int32(0))
0023 ,TimeModule = cms.untracked.PSet(limit = cms.untracked.int32(0))
0024 ,TimeEvent = cms.untracked.PSet(limit = cms.untracked.int32(0))
0025 ,MemoryCheck = cms.untracked.PSet(limit = cms.untracked.int32(0))
0026 #TimeModule, TimeEvent, TimeReport are written to LogAsbolute instead of LogInfo with a category
0027 #so they cannot be eliminated from any destination (!) unless one uses the summaryOnly option
0028 #in the Timing Service... at the price of silencing the output needed for the TimingReport profiling
0029 #
0030 #Then add the wanted ones:
0031 ,PhysicsList = cms.untracked.PSet(limit = cms.untracked.int32(-1))
0032 ,G4cout = cms.untracked.PSet(limit = cms.untracked.int32(-1))
0033 ,G4cerr = cms.untracked.PSet(limit = cms.untracked.int32(-1))
0034 ,SimG4CoreApplication = cms.untracked.PSet(limit = cms.untracked.int32(-1))
0035 )
0036 )
0037 3-Run any cmsDriver.py commands that entail simulation, e.g.(in CMSSW_3_1_0_pre4):
0038 cmsDriver.py MinBias.cfi -n 1 --step GEN,SIM --customise=Validation/Performance/TimeMemoryG4Info.py --eventcontent FEVTDEBUG --conditions FrontierConditions_GlobalTag,IDEAL_30X::All > & ! MinBias.log &
0039 The resulting file G4msg.log contains the dump of the G4 Particle Table. We run on it, extract the information we are interested in and we store it in 2 dictionaries:
0040 G4ParticleTable and G4ParticleTablePDG, the first one using Particle Name as keys, the second one using PDG ID as keys.
0041 The script also reads the HepPDT particle table contained in /src/SimGeneral/HepPDTESSource/data/pythiaparticle.tbl, it extracts all the information and stores it in 2 equivalent dictionaries,HepPdtTable and HepPdtTablePDG.
0042 Then some comparisons are made.
0043 Starting from the HepPdtTablePDG (since in CMSSW the handshake is via the PDG ID code) and look if that PDG ID exists in the G4. If it does we then check:
0044 1-For particles with ctau>10mm we check if G4 has decay tables for them (as it should)
0045 2-For all particles (that matched PDG ID) we check that they have:
0046 a-the same mass
0047 b-the same charge
0048 c-the same ctau
0049 3-For the ones that don''t match we dump some output files with the information
0050 We also dump the list of particles in HepPdtTablePDG that do not find a PDG ID match in G4ParticleTablePDG.
0051
0052 '''
0053 G4msg=open('G4msg.log','r')
0054
0055
0056 def WriteOut(myfile,mylist):
0057
0058 mylist=map(lambda a:str(a)+' ',mylist)
0059 mylist+='\n'
0060 myfile.writelines(mylist)
0061
0062 def WriteOutHtml(myfile,mylist):
0063
0064 myoutputline=['<tr>']
0065 mylist=map(lambda a:'<td align="center">'+str(a)+' '+'</td>',mylist)
0066 mylist+=['<tr>\n']
0067 myoutputline+=mylist
0068 myfile.writelines(mylist)
0069
0070 G4cout=False
0071 NewParticle=False
0072 G4ParticleTable={}
0073 G4ParticleTablePDG={}
0074 ParticleName=''
0075 ParticleInfoDict={}
0076 pdgcode=0
0077
0078 PDGcode=''
0079
0080 c=299.792458
0081 Tokens=[]
0082 for record in G4msg:
0083 if G4cout:
0084 if record[0:4]!='%MSG':
0085 if '--- G4ParticleDefinition ---' in record:
0086 NewParticle=True
0087 if NewParticle:
0088
0089 G4ParticleTable.update({ParticleName:ParticleInfoDict})
0090
0091 G4ParticleTablePDG.update({pdgcode:ParticleInfoDict})
0092
0093
0094 ParticleInfoDict={}
0095 NewParticle=False
0096 tokens=map(lambda a:a.strip(' \n'),record.split(':'))
0097 if 'Particle Name' in tokens[0]:
0098 ParticleName=tokens[1]
0099 ParticleInfoDict.update({'Particle Name':ParticleName})
0100 if 'PDG particle code' in tokens[0]:
0101 pdgcode=int(tokens[1].split()[0])
0102 ParticleInfoDict.update({'PDG ID':pdgcode})
0103 if 'Mass [GeV/c2]' in tokens[0]:
0104 mass=float(tokens[1].split()[0])
0105 ParticleInfoDict.update({'Mass [GeV/c2]':mass})
0106 if 'Charge [e]' in tokens[0]:
0107 charge=float(tokens[1])
0108 ParticleInfoDict.update({'Charge [e]':charge})
0109 if 'Lifetime [nsec]' in tokens[0]:
0110 lifetime=float(tokens[1])
0111 ParticleInfoDict.update({'Lifetime [nsec]':lifetime})
0112 if lifetime!=-1:
0113 ctau=c*lifetime
0114 ParticleInfoDict.update({'ctau [mm]':ctau})
0115 else:
0116 ParticleInfoDict.update({'ctau [mm]':-999999999})
0117 if 'G4DecayTable' in tokens[0]:
0118 decaytable=True
0119 ParticleInfoDict.update({'G4 Decay Table':decaytable})
0120 if 'Decay Table is not defined !!' in tokens[0]:
0121 decaytable=False
0122 ParticleInfoDict.update({'G4 Decay Table':decaytable})
0123 elif record[0:5]=='%MSG ':
0124
0125 G4cout=False
0126 if record[0:13]=='%MSG-i G4cout':
0127
0128 G4cout=True
0129 import os
0130 import math
0131 HepPdtTable={}
0132 HepPdtTablePDG={}
0133 for record in open(os.environ.data['CMSSW_RELEASE_BASE']+'/src/SimGeneral/HepPDTESSource/data/pythiaparticle.tbl','r'):
0134 if '//' in record:
0135 pass
0136 else:
0137 tokens = record.split()
0138
0139 if len(tokens)==6:
0140
0141 HepPdtTable.update({tokens[1]:{'Particle Name':tokens[1],'PDG ID':int(tokens[0]),'Charge [e]': float(tokens[2])/3, 'Mass [GeV/c2]':float(tokens[3]), 'ctau [mm]':float(tokens[5])}})
0142 HepPdtTablePDG.update({int(tokens[0]):{'Particle Name':tokens[1],'PDG ID':int(tokens[0]),'Charge [e]': float(tokens[2])/3, 'Mass [GeV/c2]':float(tokens[3]), 'ctau [mm]':float(tokens[5])}})
0143
0144
0145
0146
0147 print('Popping the first empty element of G4ParticleTable: %s '%G4ParticleTable.pop(''))
0148 print('G4ParticleTable contains ',len(G4ParticleTable),'elements')
0149 print('G4ParticleTablePDG contains ',len(G4ParticleTablePDG),'elements')
0150 if len(G4ParticleTable)>len(G4ParticleTablePDG):
0151 print("The following values were in the G4ParticleTable dictionary but not in the G4ParticleTablePDG one (multiple entries with different names but the same PDG code):")
0152 for value in G4ParticleTable.values():
0153 if value not in G4ParticleTablePDG.values():
0154 print(value)
0155 elif len(G4ParticleTablePDG)>len(G4ParticleTable):
0156 print("The following values were in the G4ParticleTablePDG dictionary but not in the G4ParticleTable one (multiple entries with different PDG codes but the same Particle Name):")
0157 for value in G4ParticleTablePDG.values():
0158 if value not in G4ParticleTable.values():
0159 print(value)
0160 print('HepPdtTable contains ',len(HepPdtTable),'elements')
0161 print('HepPdtTablePDG contains ',len(HepPdtTablePDG),'elements')
0162 if len(HepPdtTable)>len(HepPdtTablePDG):
0163 print("The following values were in the HepPdtTable dictionary but not in the HepPdtTablePDG one (multiple entries with different names but the same PDG code):")
0164 for value in HepPdtTable.values():
0165 if value not in HepPdtTablePDG.values():
0166 print(value)
0167 elif len(HepPdtTablePDG)>len(HepPdtTable):
0168 print("The following values were in the HepPdtTablePDG dictionary but not in the HepPdtTable one (multiple entries with different PDG codes but the same Particle Name):")
0169 for value in HepPdtTablePDG.values():
0170 if value not in HepPdtTable.values():
0171 print(value)
0172
0173
0174
0175
0176
0177
0178
0179 MatchingPDG=[]
0180 MatchingPDGMass=[]
0181 MassDiffPercentList=[]
0182 ctauDiffPercentList=[]
0183 NotMatchingPDG=[]
0184
0185 MatchingPDGDecayFileHtml=open('MatchingPDGDecay.html','w')
0186 MatchingPDGDecayFileHtml.writelines(['<html>\n','<body>\n','<table align="center", border=2>\n'])
0187 WriteOutHtml(MatchingPDGDecayFileHtml,['PDG ID', 'HepPdt Particle Name','G4 Particle Name',' HepPdt ctau [mm]','G4 ctau [mm]',' HepPdt Mass [GeV/c2]','G4 Mass [GeV/c2]',' G4 Decay Table'])
0188 MatchingPDGChargeFileHtml=open('MatchingPDGCharge.html','w')
0189 MatchingPDGChargeFileHtml.writelines(['<html>\n','<body>\n','<table align="center", border=2>\n'])
0190 WriteOutHtml(MatchingPDGChargeFileHtml,['PDG ID', 'HepPdt Particle Name','G4 Particle Name',' HepPdt Charge [e]','G4 Charge [e]'])
0191 NoLifeTimePDGFileHtml=open('NoLifeTimePDGG4.html','w')
0192 NoLifeTimePDGFileHtml.writelines(['<html>\n','<body>\n','<table align="center", border=2>\n'])
0193 WriteOutHtml(NoLifeTimePDGFileHtml,['PDG ID', 'HepPdt Particle Name','G4 Particle Name','HepPdt ctau [mm]','G4 ctau [mm]',' HepPdt Mass [GeV/c2]','G4 Mass [GeV/c2]'])
0194
0195
0196
0197
0198
0199
0200
0201
0202 for pdgcode in HepPdtTablePDG.keys():
0203 if pdgcode in G4ParticleTablePDG.keys():
0204
0205 MatchingPDG+=[HepPdtTablePDG[pdgcode]['Particle Name']]
0206
0207
0208
0209
0210 if G4ParticleTablePDG[pdgcode]['ctau [mm]']>10.0 or HepPdtTablePDG[pdgcode]['ctau [mm]']>10.0:
0211 if not G4ParticleTablePDG[pdgcode]['G4 Decay Table']:
0212 print("****Uh Oh No G4 Decay Table for ", G4ParticleTablePDG[pdgcode]['Particle Name'])
0213 else:
0214 WriteOutHtml(MatchingPDGDecayFileHtml,[pdgcode,HepPdtTablePDG[pdgcode]['Particle Name'],G4ParticleTablePDG[pdgcode]['Particle Name'],HepPdtTablePDG[pdgcode]['ctau [mm]'],G4ParticleTablePDG[pdgcode]['ctau [mm]'],HepPdtTablePDG[pdgcode]['Mass [GeV/c2]'],G4ParticleTablePDG[pdgcode]['Mass [GeV/c2]'],G4ParticleTablePDG[pdgcode]['G4 Decay Table']])
0215
0216 if G4ParticleTablePDG[pdgcode]['ctau [mm]']==-999999999:
0217 WriteOutHtml(NoLifeTimePDGFileHtml,[pdgcode,HepPdtTablePDG[pdgcode]['Particle Name'],G4ParticleTablePDG[pdgcode]['Particle Name'],HepPdtTablePDG[pdgcode]['ctau [mm]'],G4ParticleTablePDG[pdgcode]['ctau [mm]'],HepPdtTablePDG[pdgcode]['Mass [GeV/c2]'],G4ParticleTablePDG[pdgcode]['Mass [GeV/c2]']])
0218 G4ParticleTablePDG[pdgcode]['ctau [mm]']=0.0
0219 if HepPdtTablePDG[pdgcode]['ctau [mm]']!=G4ParticleTablePDG[pdgcode]['ctau [mm]']:
0220 ctauDiff=HepPdtTablePDG[pdgcode]['ctau [mm]']-G4ParticleTablePDG[pdgcode]['ctau [mm]']
0221 if HepPdtTablePDG[pdgcode]['ctau [mm]']!=0:
0222 ctauDiffPercent=math.fabs(ctauDiff/HepPdtTablePDG[pdgcode]['ctau [mm]']*100)
0223 ctauDiffPercentList+=[(abs(ctauDiffPercent),pdgcode,ctauDiff,ctauDiffPercent)]
0224 elif G4ParticleTablePDG[pdgcode]['ctau [mm]']!=0:
0225 ctauDiffPercent=math.fabs(ctauDiff/G4ParticleTablePDG[pdgcode]['ctau [mm]']*100)
0226 ctauDiffPercentList+=[(abs(ctauDiffPercent),pdgcode,ctauDiff,ctauDiffPercent)]
0227 else:
0228 ctauDiffPercent=0.0
0229 ctauDiffPercentList+=[(abs(ctauDiffPercent),pdgcode,ctauDiff,ctauDiffPercent)]
0230
0231
0232
0233 if HepPdtTablePDG[pdgcode]['Mass [GeV/c2]']!=G4ParticleTablePDG[pdgcode]['Mass [GeV/c2]']:
0234 MassDiff=HepPdtTablePDG[pdgcode]['Mass [GeV/c2]']-G4ParticleTablePDG[pdgcode]['Mass [GeV/c2]']
0235 MassDiffPercent=math.fabs(MassDiff/HepPdtTablePDG[pdgcode]['Mass [GeV/c2]']*100)
0236 MassDiffPercentList+=[(abs(MassDiffPercent),pdgcode,MassDiff,MassDiffPercent)]
0237 print(pdgcode,HepPdtTablePDG[pdgcode]['Particle Name'],G4ParticleTablePDG[pdgcode]['Particle Name'],' Mass:',HepPdtTablePDG[pdgcode]['Mass [GeV/c2]'],' Mass:',G4ParticleTablePDG[pdgcode]['Mass [GeV/c2]'],MassDiff,MassDiffPercent)
0238
0239 else:
0240
0241 MatchingPDGMass+=[HepPdtTablePDG[pdgcode]['Particle Name']]
0242
0243 if HepPdtTablePDG[pdgcode]['Charge [e]']!=G4ParticleTablePDG[pdgcode]['Charge [e]']:
0244 print(pdgcode,HepPdtTablePDG[pdgcode]['Particle Name'],G4ParticleTablePDG[pdgcode]['Particle Name'],' Charge [e]:',HepPdtTablePDG[pdgcode]['Charge [e]'],' Charge [e]:',G4ParticleTablePDG[pdgcode]['Charge [e]'])
0245 WriteOutHtml(MatchingPDGChargeFileHtml,[pdgcode,HepPdtTablePDG[pdgcode]['Particle Name'],G4ParticleTablePDG[pdgcode]['Particle Name'],HepPdtTablePDG[pdgcode]['Charge [e]'],G4ParticleTablePDG[pdgcode]['Charge [e]']])
0246 else:
0247
0248
0249
0250 NotMatchingPDG+=[(HepPdtTablePDG[pdgcode]['ctau [mm]'],pdgcode,HepPdtTablePDG[pdgcode]['Particle Name'])]
0251
0252
0253
0254
0255 MatchingPDGMassSortedFileHtml=open('MatchingPDGMassSorted.html','w')
0256 MatchingPDGMassSortedFileHtml.writelines(['<html>\n','<body>\n','<table align="center", border=2>\n'])
0257 WriteOutHtml(MatchingPDGMassSortedFileHtml,['PDG ID', 'HepPdt Particle Name','G4 Particle Name',' HepPdt Mass [GeV/c2]','G4 Mass [GeV/c2]','Mass Diff [GeV/c2]','Mass Diff %','G4 ctau [mm]'])
0258 MassDiffPercentList.sort(reverse=True)
0259 for element in MassDiffPercentList:
0260 pdgcode=element[1]
0261 MassDiff=element[2]
0262 MassDiffPercent=element[3]
0263 WriteOutHtml(MatchingPDGMassSortedFileHtml,[pdgcode,HepPdtTablePDG[pdgcode]['Particle Name'],G4ParticleTablePDG[pdgcode]['Particle Name'],HepPdtTablePDG[pdgcode]['Mass [GeV/c2]'],G4ParticleTablePDG[pdgcode]['Mass [GeV/c2]'],MassDiff,MassDiffPercent,G4ParticleTablePDG[pdgcode]['ctau [mm]']])
0264
0265
0266 MatchingPDGctauSortedFileHtml=open('MatchingPDGctauSorted.html','w')
0267 MatchingPDGctauSortedFileHtml.writelines(['<html>\n','<body>\n','<table align="center", border=2>\n'])
0268 WriteOutHtml(MatchingPDGctauSortedFileHtml,['PDG ID', 'HepPdt Particle Name','G4 Particle Name',' HepPdt ctau [mm]','G4 ctau [mm]','ctau Diff [mm]','ctau Diff %'])
0269 ctauDiffPercentList.sort(reverse=True)
0270 for element in ctauDiffPercentList:
0271 pdgcode=element[1]
0272 ctauDiff=element[2]
0273 ctauDiffPercent=element[3]
0274 WriteOutHtml(MatchingPDGctauSortedFileHtml,[pdgcode, HepPdtTablePDG[pdgcode]['Particle Name'],G4ParticleTablePDG[pdgcode]['Particle Name'], HepPdtTablePDG[pdgcode]['ctau [mm]'],G4ParticleTablePDG[pdgcode]['ctau [mm]'],ctauDiff,ctauDiffPercent])
0275
0276
0277
0278 NotMatchingPDGFileHtml=open('NotMatchingPDG.html','w')
0279 NotMatchingPDGFileHtml.writelines(['<html>\n','<body>\n','<table align="center", border=2>\n'])
0280 WriteOutHtml(NotMatchingPDGFileHtml,['PDG ID', 'HepPdt Particle Name','HepPdt Mass [GeV/c2]',' HepPdt Charge [e]','HepPdt ctau [mm]'])
0281 NotMatchingPDG.sort(reverse=True)
0282 for element in NotMatchingPDG:
0283 pdgcode=element[1]
0284 WriteOutHtml(NotMatchingPDGFileHtml,[pdgcode,HepPdtTablePDG[pdgcode]['Particle Name'],HepPdtTablePDG[pdgcode]['Mass [GeV/c2]'],HepPdtTablePDG[pdgcode]['Charge [e]'],HepPdtTablePDG[pdgcode]['ctau [mm]']])
0285
0286
0287 MatchingPDGDecayFileHtml.writelines(['</table>\n','<body>\n','</html>\n'])
0288 MatchingPDGChargeFileHtml.writelines(['</table>\n','<body>\n','</html>\n'])
0289 MatchingPDGMassSortedFileHtml.writelines(['</table>\n','<body>\n','</html>\n'])
0290 MatchingPDGctauSortedFileHtml.writelines(['</table>\n','<body>\n','</html>\n'])
0291 NotMatchingPDGFileHtml.writelines(['</table>\n','<body>\n','</html>\n'])
0292 NoLifeTimePDGFileHtml.writelines(['</table>\n','<body>\n','</html>\n'])
0293
0294
0295