iceprod.modules.dcorsika

35 """ 36 This class provides an interface for preprocessing files in iceprod 37 """ 38

39 - def __init__(self):

40 IPBaseClass.__init__(self) 41 self.logger = logging.getLogger('iceprod::Corsika') 42 self.name = 'corsika' 43 self.parameters['arrang'] = -119. 44 45 self.AddParameter('version','Corsika version','v6720') 46 self.AddParameter('platform','compliler platform','') 47 self.AddParameter('cache','Should cache taball?',False) 48 self.AddParameter('cachedir','Cache directory', 49 '$system(cache)/%s_icesoft/%s'% (getpass.getuser(), self.name)) 50 self.AddParameter('URL','fetch tarball from URL',None) 51 self.AddParameter('runnum','Run number','') 52 self.AddParameter('seed','Random seed','1') 53 self.AddParameter('nevents','Number of Events',0) 54 self.AddParameter('outfile','Output file',"DAT%(runnum)06d") 55 self.AddParameter('logfile','Log file','%(outfile)s.log') 56 self.AddParameter('outdir','Output directory','.') 57 self.AddParameter('topdir','Top directory','') 58 self.AddParameter('tmpdir','Temporary directory','%(outdir)s/dcors%(runnum)d') 59 self.AddParameter('model','Physics Model','SIBYLL') 60 self.AddParameter('lemodel','Low Energy Physics Model','gheisha') 61 self.AddParameter('donkg','Run NKG',0) # no NKG by default 62 self.AddParameter('doegs','Run EGS',0) # no EGS by default 63 self.AddParameter('eslope','CR spectral index (only if ranpri=0)',-2.7) 64 self.AddParameter('crtype','CR primary type',14) 65 self.AddParameter('cthmin','Min theta of injected cosmic rays',0.0) 66 self.AddParameter('cthmax','Max theta of injected cosmic rays',89.99) 67 self.AddParameter('cphmin','Min phi of injected cosmic rays',0.0) 68 self.AddParameter('cphmax','Max phi of injected cosmic rays',360.0) 69 self.AddParameter('emin','CR min energy',600.) 70 self.AddParameter('emax','CR max energy',1.e11) 71 self.AddParameter('ecuts','hadron/em energy cut (deprecated: instead use ecuts(i),i=1..4)',0) 72 self.AddParameter('ecuts1','hadron min energy (see corsika docs)',273) 73 self.AddParameter('ecuts2','muon min energy (see corsika docs)',273) 74 self.AddParameter('ecuts3','electron min energy (see corsika docs)',0.003) 75 self.AddParameter('ecuts4','photon min energy (see corsika docs)',0.003) 76 self.AddParameter('atmod','Atmosphere model (October=13)',13) 77 self.AddParameter('debug','boolean: enable disable debug mode',False)

78 79

80 - def stage(self):

81 """ 82 Stage files and executables 83 """ 84 85 par = self.parameters 86 par['tmpdir'] = self.GetParameter('tmpdir') % par; 87 88 self.logger.info('setting up working directory: %(tmpdir)s' % par) 89 if not os.path.exists(par['tmpdir']): 90 os.makedirs(par['tmpdir']); # create temporary directory 91 os.chdir(par['tmpdir']); # cd to temporary directory 92 93 self.logger.info('caching: %s' % self.GetParameter('cache')) 94 if self.GetParameter('cache'): 95 cachedir = self.parser.parse(self.GetParameter('cachedir')) 96 if not os.path.exists(cachedir): 97 os.makedirs(cachedir) 98 else: 99 cachedir = expandvars("$PWD") 100 baseurl = par['url'] 101 102 meta = I3Tarball() 103 meta.version = self.GetParameter('version') 104 meta.platform = self.GetParameter('platform') 105 meta.name = self.name 106 meta.suffix = ".tar.gz" 107 if meta.platform: 108 meta.filebase = "%s-%s.%s" % (meta.name, meta.version, meta.platform) 109 else: 110 meta.filebase = "%s-%s" % (meta.name, meta.version) 111 meta.md5sum = "%s.md5sum" % meta.filebase 112 meta.url = "%s/%s" % (baseurl,meta.filebase) 113 exe.fetch_tarball(meta,cachedir) 114 par['topdir'] = meta.path 115 116 # link data files for corsika 117 # Files necessary for QGSJET and QGSJET-II included 118 # DPMJET and VENUS files are *not* 119 os.symlink("%(topdir)s/bin/NUCNUCCS" % par, "NUCNUCCS"); 120 os.symlink("%(topdir)s/bin/QGSDAT01" % par, "QGSDAT01"); 121 os.symlink("%(topdir)s/bin/SECTNU" % par, "SECTNU"); 122 os.symlink("%(topdir)s/bin/qgsdat-II-03" % par, "qgsdat-II-03"); 123 os.symlink("%(topdir)s/bin/sectnu-II-03" % par, "sectnu-II-03"); 124 os.symlink("%(topdir)s/bin/GLAUBTAR.DAT" % par, "GLAUBTAR.DAT"); 125 os.symlink("%(topdir)s/bin/NUCLEAR.BIN" % par, "NUCLEAR.BIN"); 126 os.environ['LD_LIBRARY_PATH'] = expandvars("%(topdir)s/lib:$LD_LIBRARY_PATH" % par); 127 128 # If we have EGS files too, link those 129 egsfiles = glob.glob("%(topdir)s/bin/EGSDAT*" % par) 130 131 # If we have EPOS files too, link those 132 eposfiles = glob.glob("%(topdir)s/bin/epos.*" % par) 133 134 for file in egsfiles + eposfiles: 135 os.symlink(file, os.path.basename(file)); 136 137 # If FLUKA exists, use it 138 if os.path.exists("%(topdir)s/fluka" % par): 139 os.environ['FLUPRO'] = expandvars("%(topdir)s/fluka" % par); 140 141 return meta.path

142 143 144

145 - def Execute(self,stats):

146 """ 147 Run CORSIKA 148 corsika output is stdout: must create a temporary directory and cd there 149 """ 150 151 cwd = os.getcwd() 152 par = self.parameters 153 154 # Retrieve tarball and stage environment 155 self.stage() 156 self.configure() 157 self.write_steering() 158 159 # CORSIKA binary 160 # New standard binary name style 161 par['model'] = par['model'].upper() 162 par['versionnumber'] = par['version'].lstrip('v') 163 par['corsbinary'] = "%(topdir)s/bin/corsika%(versionnumber)sLinux_%(model)s_%(lemodel)s" % par; 164 if not os.path.exists(par['corsbinary']): 165 os.chdir(cwd); # cd back to original directory 166 print >> sys.stderr,"CORSIKA binary does not exist: corsika%(versionnumber)sLinux_%(model)s_%(lemodel)s\n" % par; 167 raise Exception,"CORSIKA binary does not exist: corsika%(versionnumber)sLinux_%(model)s_%(lemodel)s\n" % par; 168 # Old symlink style 169 #par['corsbinary'] = "%(topdir)s/bin/corsika.%(model)s.Linux" % par; 170 system("cp %(corsbinary)s corsika.%(runnum)s.Linux" % par); 171 172 # Corsika output file 173 par['outfile'] = os.path.join(expandvars(par['outdir']%par),expandvars(par['outfile']%par)) 174 par['corout'] = "DAT%(runnum)06d" % par; # Plain CORSIKA output 175 par['logfile'] = par['logfile'] % par 176 177 # Execution command 178 cors_cmd = "%(tmpdir)s/corsika.%(runnum)d.Linux < %(inputfile)s > %(logfile)s " % par; 179 self.logger.info(cors_cmd); 180 181 # Run CORSIKA 182 status, output = getstatusoutput(cors_cmd); 183 self.logger.info(output); 184 if status: 185 os.chdir(cwd); # cd back to original directory 186 raise Exception, "dCorsika python sucks! %s\n" % output; 187 188 # Check if dCORSIKA Output exists 189 if os.path.exists("%(tmpdir)s/%(corout)s.gz" % par): 190 system("cp %(tmpdir)s/%(corout)s.gz %(outfile)s" % par); 191 # Check if CORSIKA Output exists 192 elif os.path.exists("%(tmpdir)s/%(corout)s" % par): 193 system("cp %(tmpdir)s/%(corout)s %(outfile)s" % par); 194 # if neither output exists, quit 195 else: 196 os.chdir(cwd); # cd back to original directory 197 print >> sys.stderr,"CORSIKA Output does not exist: exit before UCR\n"; 198 raise Exception,"CORSIKA Output does not exist: exit before UCR\n"; 199 200 201 202 # check if the output is OK 203 nevcorsika = 0 204 try: 205 status,nevcorsika=getstatusoutput("cat %(logfile)s|grep \"GENERATED EVENTS\"|grep -v NEV|awk '{ print $6 }'" % par); 206 nevcorsika = int(nevcorsika.strip()); 207 stats['nevcorsika'] = nevcorsika 208 except Exception,e: 209 self.logger.error(e); 210 if nevcorsika == self.GetParameter('nevents'): 211 system('touch %(outdir)s/corsika.%(runnum)d.isok' % par) 212 self.logger.info("OK"); 213 print >> sys.stderr, "Corsika OK\n"; 214 else : 215 system('touch %(outdir)s/corsika.%(runnum)d.isnotok' % par) 216 print >> sys.stderr, "Corsika not OK\n"; 217 self.logger.error("NOT OK"); 218 self.logger.error(exe.tail(self.GetParameter('logfile'))); 219 os.chdir(cwd); # cd back to original directory 220 return 1; 221 222 os.chdir(cwd); # cd back to original directory 223 return 0;

224

225 - def configure(self):

226 """ 227 Configure and write INPUTS steering file 228 """ 229 par = self.parameters 230 seed = self.GetParameter('seed') 231 par['seed1'] = int(seed)+0; 232 par['seed2'] = int(seed)+1; 233 par['seed3'] = int(seed)+2; 234 235 # NKG/EGS 236 NKGparams = ""; 237 nkg = 'F' 238 egs = 'F' 239 if par['donkg']: nkg = 'T' 240 if par['doegs']: egs = 'T' 241 NKGparams = "ELMFLG %s %s " % (nkg,egs) 242 NKGparams += " em. interaction flags (NKG,EGS)\n"; 243 if par['donkg']: # if NKG parameterizations ON 244 NKGparams += "RADNKG 2.E5 outer radius for NKG lat.dens.determ."; 245 par['NKGparams'] = NKGparams 246 247 # Construct HE interaction model steering commands 248 modelStr = ""; 249 model = self.GetParameter('model') 250 if model in ("qgsjet","qgsii"): 251 modelStr += "QGSJET T 0 use qgsjet for high energy hadrons\n"; 252 modelStr += "QGSSIG T use qgsjet hadronic cross sections"; 253 elif model == "dpmjet": 254 modelStr += "DPMJET T 0 use dpmjet for high energy hadrons\n"; 255 modelStr += "DPJSIG T all hail Glaubtar!"; 256 elif model == "sibyll": 257 modelStr += "SIBYLL T 0 use sibyll for high energy hadrons\n"; 258 modelStr += "SIBSIG T use sibyll hadronic cross sections"; 259 elif model == "epos": 260 modelStr += "EPOS T 0 use epos for high energy hadrons\n"; 261 modelStr += "EPOSIG T use epos hadronic cross sections\n"; 262 modelStr += "EPOPAR input epos.param !initialization input file for epos\n"; 263 modelStr += "EPOPAR fname inics epos.inics !initialization input file for epos\n"; 264 modelStr += "EPOPAR fname iniev epos.iniev !initialization input file for epos\n"; 265 modelStr += "EPOPAR fname initl epos.initl !initialization input file for epos\n"; 266 modelStr += "EPOPAR fname inirj epos.inirj !initialization input file for epos\n"; 267 modelStr += "EPOPAR fname inihy epos.ini1b !initialization input file for epos\n"; 268 modelStr += "EPOPAR fname check none !dummy output file for epos\n"; 269 modelStr += "EPOPAR fname histo none !dummy output file for epos\n"; 270 modelStr += "EPOPAR fname data none !dummy output file for epos\n"; 271 modelStr += "EPOPAR fname copy none !dummy output file for epos"; 272 273 # Turn on/off dCORSIKA debugging 274 if self.GetParameter('debug'): 275 par['debug'] = "T"; 276 else: 277 par['debug'] = "F"; 278 279 # Check if old-style ecuts parameter is set 280 ecuts = self.GetParameter('ecuts') 281 if ecuts: 282 self.SetParameter('ecuts1',ecuts) 283 self.SetParameter('ecuts2',ecuts) 284 285 # Convert input phi from IceCube Coordinates to CORSIKA Coordinates 286 # CORSIKA will rotate the particles back to IceCube Coordinates in 287 # the output routine. 288 par['cphmin_cc'] = par['cphmin'] + par['arrang'] 289 par['cphmax_cc'] = par['cphmax'] + par['arrang'] 290 291 # Check the domain of min and max phi and fix if needed 292 if par['cphmax_cc'] - par['cphmin_cc'] > 360.0: 293 self.logger.error('Phi range greater than 360deg') 294 if par['cphmin_cc'] < -360.0: 295 par['cphmin_cc'] += 360.0 296 par['cphmax_cc'] += 360.0 297 if par['cphmax_cc'] > 360.0: 298 par['cphmin_cc'] -= 360.0 299 par['cphmax_cc'] -= 360.0 300 301 302 # Write out the Corsika INPUTS steering file 303 input = "" 304 input += "RUNNR %(runnum)d number of run\n"; 305 input += "EVTNR 1 number of first shower event\n"; 306 input += "NSHOW %(nevents)d number of showers to generate\n"; 307 input += "PRMPAR %(crtype)d particle type of prim. particle\n"; 308 input += "ESLOPE %(eslope)f slope of primary energy spectrum\n"; 309 input += "ERANGE %(emin)f %(emax)f energy range of primary particle\n"; 310 input += "THETAP %(cthmin)f %(cthmax)f range of zenith angle (degree)\n"; 311 input += "PHIP %(cphmin_cc)f %(cphmax_cc)f range of azimuth angle (degree)\n"; 312 input += "SEED %(seed1)d 0 0 seed for 1. random number sequence\n"; 313 input += "SEED %(seed2)d 0 0 seed for 2. random number sequence\n"; 314 input += "SEED %(seed3)d 0 0 seed for 3. random number sequence\n"; 315 input += "OBSLEV 2834.E2 observation level (in cm)\n"; 316 input += "%(NKGparams)s \n"; 317 input += "ARRANG %(arrang)f rotation of array to north\n"; 318 input += "FIXHEI 0. 0 first interaction height & target\n"; 319 input += "FIXCHI 0. starting altitude (g/cm**2)\n"; 320 input += "MAGNET 16.4 -53.4 magnetic field south pole\n"; 321 input += "HADFLG 0 1 0 1 0 2 flags hadr.interact. & fragmentation\n"; 322 input += "%s \n" % modelStr; 323 input += "ECUTS %(ecuts1).04f %(ecuts2).04f %(ecuts3).04f %(ecuts4).04f energy cuts for particles\n"; 324 input += "MUADDI T additional info for muons\n"; 325 input += "MUMULT T muon multiple scattering angle\n"; 326 input += "LONGI F 20. F F longit.distr. & step size & fit\n"; 327 input += "MAXPRT 0 max. number of printed events\n"; 328 input += "ECTMAP 100 cut on gamma factor for printout\n"; 329 input += "STEPFC 1.0 mult. scattering step length fact.\n"; 330 input += "DEBUG %(debug)s 6 F 1000000 debug flag and log.unit for out\n"; 331 input += "DIRECT ./ output directory\n"; 332 input += "ATMOD %(atmod)s october atmosphere\n"; 333 self.steering = input

334 335

336 - def write_steering(self):

337 # terminate steering and write it to file 338 par = self.parameters 339 self.steering += "EXIT terminates input\n"; 340 341 tmpdir = par['tmpdir'] 342 inputfile = open(tmpdir+"/INPUTS",'w'); 343 inputfile.write(self.steering % par); 344 par['inputfile'] = inputfile.name; 345 inputfile.close();

Source Code for Module iceprod.modules.dcorsika