| 110 | | def do_add_time_of_flight(self, line): |
| 111 | | print("Running patched do_add_time_of_flight ") |
| 112 | | |
| 113 | | args = self.split_arg(line) |
| 114 | | # check the validity of the arguments and reformat args |
| 115 | | self. check_add_time_of_flight(args) |
| 116 | | |
| 117 | | event_path, threshold = args |
| 118 | | # gunzip the file |
| 119 | | if event_path.endswith('.gz'): |
| 120 | | need_zip = True |
| 121 | | misc.gunzip(event_path) |
| 122 | | event_path = event_path[:-3] |
| 123 | | else: |
| 124 | | need_zip = False |
| 125 | | |
| 126 | | import random |
| 127 | | try: |
| 128 | | import madgraph.various.lhe_parser as lhe_parser |
| 129 | | except: |
| 130 | | import internal.lhe_parser as lhe_parser |
| 131 | | |
| 132 | | logger.info('Add time of flight information on file %s' % event_path) |
| 133 | | lhe = lhe_parser.EventFile(event_path) |
| 134 | | output = open('%s_2vertex.lhe' % event_path, 'w') |
| 135 | | # write the banner to the output file |
| 136 | | output.write(lhe.banner) |
| 137 | | |
| 138 | | # get the associate param_card |
| 139 | | begin_param = lhe.banner.find('<slha >') |
| 140 | | end_param = lhe.banner.find('</slha >') |
| 141 | | param_card = lhe.banner[begin_param + 6: end_param]. split('\n') |
| 142 | | param_card = param_card_mod.ParamCard(param_card) |
| 143 | | |
| 144 | | cst = 6.58211915e-25 # hbar in GeV s |
| 145 | | c = 299792458000 # speed of light in mm/s |
| 146 | | sm_lepton_list = [11, 12, 13, 14, 15, 16] |
| 147 | | |
| 148 | | pspss_n_list = [8000011, 8000012] |
| 149 | | mass_splitting = param_card.get_value('PSPSS ', 2) |
| 150 | | damping = param_card.get_value('PSPSS ', 6) |
| 151 | | |
| 152 | | # Loop over all events |
| 153 | | for event in lhe: |
| 154 | | leptonnumber = 0 |
| 155 | | for particle in event: |
| 156 | | if particle.status == 1: |
| 157 | | if particle.pid in sm_lepton_list: |
| 158 | | leptonnumber += 1 |
| 159 | | elif -particle.pid in sm_lepton_list: |
| 160 | | leptonnumber -= 1 |
| 161 | | |
| 162 | | write_event = True |
| 163 | | for particle in event: |
| 164 | | id = particle.pid |
| 165 | | width = param_card['decay '].get((abs(id),)).value |
| 166 | | |
| 167 | | if width: |
| 168 | | if id in pspss_n_list: |
| 169 | | tau0 = random.expovariate(width / cst) |
| 170 | | |
| 171 | | if 0.5 * (1 + math.exp(-damping) * math.cos(mass_splitting * tau0 / cst)) >= random.random(): |
| 172 | | write_event = (leptonnumber == 0) |
| 173 | | else: |
| 174 | | write_event = (leptonnumber != 0) |
| 175 | | |
| 176 | | else: |
| 177 | | tau0 = random.expovariate(width / cst) |
| 178 | | |
| 179 | | vtim = c * tau0 |
| 180 | | if vtim > threshold: |
| 181 | | particle.vtim = vtim |
| 182 | | |
| 183 | | # write this modify event |
| 184 | | if write_event: |
| 185 | | output.write(str(event)) |
| 186 | | |
| 187 | | |
| 188 | | output.write('</LesHouchesEvents >\n') |
| 189 | | output.close() |
| 190 | | files.mv('%s_2vertex.lhe' % event_path, event_path) |
| 191 | | |
| 192 | | ifneed_zip: misc.gzip(event_path) |
| 193 | | |
| | 110 | def do_add_time_of_flight(self, line): |
| | 111 | print("Running patched do_add_time_of_flight") |
| | 112 | |
| | 113 | args = self.split_arg(line) |
| | 114 | #check the validity of the arguments and reformat args |
| | 115 | self.check_add_time_of_flight(args) |
| | 116 | |
| | 117 | event_path, threshold = args |
| | 118 | #gunzip the file |
| | 119 | if event_path.endswith('.gz'): |
| | 120 | need_zip = True |
| | 121 | misc.gunzip(event_path) |
| | 122 | event_path = event_path[:-3] |
| | 123 | else: |
| | 124 | need_zip = False |
| | 125 | |
| | 126 | import random |
| | 127 | try: |
| | 128 | import madgraph.various.lhe_parser as lhe_parser |
| | 129 | except: |
| | 130 | import internal.lhe_parser as lhe_parser |
| | 131 | |
| | 132 | logger.info('Add time of flight information on file %s' % event_path) |
| | 133 | lhe = lhe_parser.EventFile(event_path) |
| | 134 | output = open('%s_2vertex.lhe' % event_path, 'w') |
| | 135 | #write the banner to the output file |
| | 136 | output.write(lhe.banner) |
| | 137 | |
| | 138 | # get the associate param_card |
| | 139 | begin_param = lhe.banner.find('<slha>') |
| | 140 | end_param = lhe.banner.find('</slha>') |
| | 141 | param_card = lhe.banner[begin_param+6:end_param].split('\n') |
| | 142 | param_card = param_card_mod.ParamCard(param_card) |
| | 143 | |
| | 144 | cst = 6.58211915e-25 # hbar in GeV s |
| | 145 | c = 299792458000 # speed of light in mm/s |
| | 146 | sm_lepton_list = [11, 12, 13, 14, 15, 16] |
| | 147 | |
| | 148 | pspss_n_list = [8000011, 8000012] |
| | 149 | mass_splitting = param_card.get_value('PSPSS', 2) |
| | 150 | damping = param_card.get_value('PSPSS', 6) |
| | 151 | |
| | 152 | # Loop over all events |
| | 153 | for event in lhe: |
| | 154 | leptonnumber = 0 |
| | 155 | for particle in event: |
| | 156 | if particle.status == 1: |
| | 157 | if particle.pid in sm_lepton_list: |
| | 158 | leptonnumber += 1 |
| | 159 | elif -particle.pid in sm_lepton_list: |
| | 160 | leptonnumber -= 1 |
| | 161 | |
| | 162 | write_event = True |
| | 163 | for particle in event: |
| | 164 | id = particle.pid |
| | 165 | width = param_card['decay'].get((abs(id),)).value |
| | 166 | |
| | 167 | if width: |
| | 168 | if id in pspss_n_list: |
| | 169 | tau0 = random.expovariate(width / cst) |
| | 170 | |
| | 171 | if 0.5 * (1 + math.exp(-damping) * math.cos(mass_splitting * tau0 / cst)) >= random.random(): |
| | 172 | write_event = (leptonnumber == 0) |
| | 173 | else: |
| | 174 | write_event = (leptonnumber != 0) |
| | 175 | |
| | 176 | else: |
| | 177 | tau0 = random.expovariate(width / cst) |
| | 178 | |
| | 179 | vtim = c * tau0 |
| | 180 | if vtim > threshold: |
| | 181 | particle.vtim = vtim |
| | 182 | |
| | 183 | # write this modify event |
| | 184 | if write_event: |
| | 185 | output.write(str(event)) |
| | 186 | |
| | 187 | output.write('</LesHouchesEvents>\n') |
| | 188 | output.close() |
| | 189 | files.mv('%s_2vertex.lhe' % event_path, event_path) |
| | 190 | |
| | 191 | if need_zip: |
| | 192 | misc.gzip(event_path) |