#!/bin/python import os, shutil from xml.dom import minidom from numpy import append, savetxt, shape, array, newaxis, zeros, arange from modules.read import ReadResults, ReadProfile, resultdir from modules.spectrum import spectrum from modules.maps import computeMap from modules.arrival_angle import arrivalAngle from modules.timing import timing from modules.observables import param_vs_energy, delay_vs_theta from modules.constants import degre xmlfile = minidom.parse("simulations.xml") simulations = xmlfile.getElementsByTagName("simu") for simu in simulations: fileId = simu.getAttribute("simulation_dir") output_dir = resultdir+simu.getAttribute("id")+"/" if not os.path.exists(output_dir): os.makedirs(output_dir) powerlaw_index = float(simu.getAttribute("powerlaw_index")) print "#=============================================================================#" print " Start working on ", fileId print " Output directory: ", output_dir print " powerlaw index: ", powerlaw_index print "#=============================================================================#" # copy profile shutil.copy(fileId+"/profile.dat",output_dir+"profile.dat") # read files print " > Reading data" n_phot, n_lept = ReadProfile("../"+fileId,[4,5]) ratio = n_phot#/n_lept results=ReadResults(fileId) # select photons only results = results[:,results[0]%2==0] print " ... ", shape(results)[1], "events selected" generation = results[0] weightini = results[1]/ratio energy = results[2] time = results[3] arrival_pos1 = results[4]*degre arrival_pos2 = results[6]*degre theta = (results[4]-results[6])*degre phi = (results[5]-results[7])*degre theta_arrival = results[8]*degre Esource =results[9] theta_min = 1e-5 theta_max = 90 max_dt = max(time) #=============================================================================# # NO SELECTION #=============================================================================# print " > No selection of events ..." # SPECTRUM (SOURCE) =========================================================# weight_source = (Esource/min(Esource))**(1-powerlaw_index) weight = weightini* weight_source nbBins = 100 Es=array(list(set(Esource))) Ws= (Es/min(Es))**(1-powerlaw_index) / ratio Es,Fs = spectrum(Es,Ws,nbBins=nbBins) Source = Es[:,newaxis] Source = append(Source,Fs[:,newaxis],axis=1) if shape(Ws!=0)[0]>1: ratio = (max(Es)-min(Es))*(max(Fs)+min(Fs))/2 weight /= ratio # IMAGING ===================================================================# print " ... Computing image" nbBins = 50 computeMap(theta,phi,weight,energy,output_dir,nbBins,borne=[theta_max,theta_max]) # SPECTRUM (MEASURED) =======================================================# print " ... Computing spectrum" nbBins = 100 ener,flux = spectrum(energy,weight,nbBins) Spectrum = ener[:,newaxis] Spectrum = append(Spectrum,flux[:,newaxis],axis=1) # ANGLE VERSUS ENERGY =======================================================# print " ... Computing arrival angle versus energy" nbBins = 100 ener,angle = param_vs_energy(theta_arrival,energy,nbBins) Angle_Energy = ener[:,newaxis] Angle_Energy = append(Angle_Energy,angle[:,newaxis],axis=1) # Arrival angle =============================================================# print " ... Computing arrival angle distribution" nbBins = 100 theta2,dndtheta = arrivalAngle(theta_arrival,weight,nbBins,theta_range=[theta_min,theta_max]) arrival_Angle = theta2[:,newaxis] arrival_Angle = append(arrival_Angle,dndtheta[:,newaxis],axis=1) # TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE =============================# print " ... Computing time distribution" nbBins = 200 delta_t,dNdt = timing(time,weight,nbBins)#,dt_range=[]) Timing = delta_t[:,newaxis] Timing = append(Timing,dNdt[:,newaxis],axis=1) angle,dt = delay_vs_theta(theta_arrival,time,nbBins)#,dt_range=[]) Delay_vs_angle = angle[:,newaxis] Delay_vs_angle = append(Delay_vs_angle,dt[:,newaxis],axis=1) #=============================================================================# # BY ENERGY RANGE #=============================================================================# print " > Selection by energy range ..." Emin = [1e-3,1e0,1e3] #GeV Emax = [1e0,1e3,1e5] #GeV for n in arange(0,3,1): cond= (energy>Emin[n]) & (energy Selection by time range ..." tmax = [2.6298e5,3.15576e7,1.57788e8] # seconds (1 month, 1 year, 5 years) for n in arange(0,3,1): cond= (time Selection by jet opening angle ..." opening_angle = [45,30,10,5,1] # degres for n in arange(0,5,1): cond= (arrival_pos1<=opening_angle[n]) & (arrival_pos2<=opening_angle[n]) # Arrival angle ==========================================================# nbBins = 100 theta2,dndtheta = arrivalAngle(theta_arrival[cond],weight[cond],nbBins,theta_range=[theta_min,theta_max]) arrival_Angle = append(arrival_Angle,dndtheta[:,newaxis],axis=1) print " ... ", ((n+1)*100)/5,"% done" #=============================================================================# # BY GENERATION #=============================================================================# print " > Selection by generation ..." NbTotEvents = sum(weight) gen_tab =list(set(generation)) Gen_cont = zeros([shape(gen_tab)[0],2]) for gen in gen_tab: cond = generation==gen contrib =sum(weight[cond])/NbTotEvents*100 i = gen_tab.index(gen) Gen_cont[i,0]=int(gen) Gen_cont[i,1]=contrib print " ... gen=",int(gen),"-> contribution:",int(contrib),"%" # SPECTRUM (MEASURED) ====================================================# nbBins = 100 ener,flux = spectrum(energy[cond],weight[cond],nbBins) Spectrum = append(Spectrum,flux[:,newaxis],axis=1) # ANGLE VERSUS ENERGY ====================================================# nbBins = 100 ener,angle = param_vs_energy(theta_arrival[cond],energy[cond],nbBins) Angle_Energy = append(Angle_Energy,angle[:,newaxis],axis=1) # arrivalAngle =================================================================# nbBins = 100 theta2,dndtheta = arrivalAngle(theta_arrival[cond],weight[cond],nbBins,theta_range=[theta_min,theta_max]) arrival_Angle = append(arrival_Angle,dndtheta[:,newaxis],axis=1) # TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE ==========================# nbBins = 200 delta_t,dNdt = timing(time[cond],weight[cond],nbBins) Timing = append(Timing,dNdt[:,newaxis],axis=1) angle,dt = delay_vs_theta(theta_arrival[cond],time[cond],nbBins)#,dt_range=[]) Delay_vs_angle = append(Delay_vs_angle,dt[:,newaxis],axis=1) #=============================================================================# print " > writing files" savetxt(output_dir+"/Spectrum.txt",Spectrum) savetxt(output_dir+"/Source_spectrum.txt",Source) savetxt(output_dir+"/Angle_vs_Energy.txt",Angle_Energy) savetxt(output_dir+"/arrival_Angle_distribution.txt",arrival_Angle) savetxt(output_dir+"/Timing.txt",Timing) savetxt(output_dir+"/Delay_versus_angle.txt",Delay_vs_angle) savetxt(output_dir+"/Generation.txt",Gen_cont) print "#=============================================================================#"