diff --git a/grand/sim/efield2voltage.py b/grand/sim/efield2voltage.py index cb9d28a7..0b30866f 100644 --- a/grand/sim/efield2voltage.py +++ b/grand/sim/efield2voltage.py @@ -314,7 +314,7 @@ def final_resample(self): self.target_lenght= int(self.target_duration_us*self.target_sampling_rate_mhz) logger.info(f"resampling the voltage from {self.f_samp_mhz[0]} to {self.target_sampling_rate_mhz} MHz, new trace lenght is {self.target_lenght} samples") #we use fourier interpolation, becouse its easy! - self.vout = sf.irfft(self.vout_f, m)*ratio*self.padding_factor #renormalize the amplitudes + self.vout = sf.irfft(self.vout_f, m)*ratio #renormalize the amplitudes #MATIAS: TODO: now, we are missing a place to store the new sampling rate! elif(self.params["add_noise"] or self.params["add_rf_chain"]): #we know we dont need to resample, but we might need to reproces the Voc (curently stored in vout by compute_voc_event) to take into acount the noise or the chain self.vout[:] = sf.irfft(self.vout_f) diff --git a/scripts/convert_efield2efield.py b/scripts/convert_efield2efield.py index efa38db2..4b1cb535 100755 --- a/scripts/convert_efield2efield.py +++ b/scripts/convert_efield2efield.py @@ -511,7 +511,7 @@ def impz(b, a): ratio=1 #to resample we use fourier interpolation, becouse it seems to be better than scipy.decimate (points are closer to the original trace) - vout = sf.irfft(vout_f, m)*ratio*padding_factor #renormalize the amplitudes. We have two effects we changed the sampling rate but we also have enlarged the trace, so the frequency content is diluted + vout = sf.irfft(vout_f, m)*ratio #renormalize the amplitudes. if(event_idx==0 and PLOT): plt.scatter(np.arange(0,len(vout[PlotAnt][0]))/ratio,vout[PlotAnt][0],label="sampled",c="red") diff --git a/sim2root/Common/IllustrateSimPipe.py b/sim2root/Common/IllustrateSimPipe.py index 3c6d1d1c..4a675e32 100644 --- a/sim2root/Common/IllustrateSimPipe.py +++ b/sim2root/Common/IllustrateSimPipe.py @@ -227,30 +227,42 @@ def plotfigure(time1="Off",signal1="Off",srate1=0,time2="Off",signal2="Off",srat #unstack the trace #I only plot the first quarter of the trace becouse that is where the peak is, and adding the long tail just puts a lot of noise on the fft, + first=600 + last=1200 - tracex=trace[du_idx,0]*10 - tracey=trace[du_idx,1]*10 - tracez=trace[du_idx,2]*10 - tracet=np.arange(0,len(tracez))*dt_ns_l0[du_idx] #+1200*dt_ns_l0[du_idx] - #print("efield",len(tracez)) - - vtracex=vtrace[du_idx,0] - vtracey=vtrace[du_idx,1] - vtracez=vtrace[du_idx,2] - vtracet=np.arange(0,len(vtracez))*dt_ns_l0[du_idx] #+1200*dt_ns_l0[du_idx] - #print("voltage",len(vtracez)) - - atracex=atrace[du_idx,0]*109.86 #this number comes from doing 0.9V/8192, the maximum of ADC divided by 13bits - atracey=atrace[du_idx,1]*109.86 - atracez=atrace[du_idx,2]*109.86 - atracet=np.arange(0,len(atracez))*dt_ns_l1[du_idx] #+300*dt_ns_l1[du_idx] - #print("adc",len(atracez)) - - rtracex=rtrace[du_idx,0]*10 - rtracey=rtrace[du_idx,1]*10 - rtracez=rtrace[du_idx,2]*10 - rtracet=np.arange(0,len(rtracez))*dt_ns_l1[du_idx] #+300*dt_ns_l1[du_idx] - #print("refield",len(rtracez)) + first_l0=int(first/dt_ns_l0[du_idx]) + first_l1=int(first/dt_ns_l1[du_idx]) + last_l0=int(last/dt_ns_l0[du_idx]) + last_l1=int(last/dt_ns_l1[du_idx]) + + + + tracex=trace[du_idx,0][first_l0:last_l0]*10 + tracey=trace[du_idx,1][first_l0:last_l0]*10 + tracez=trace[du_idx,2][first_l0:last_l0]*10 + tracet=np.arange(0,len(tracez))*dt_ns_l0[du_idx]+first + print("efield",len(tracez)) + + vtracex=vtrace[du_idx,0][first_l0:last_l0] + vtracey=vtrace[du_idx,1][first_l0:last_l0] + vtracez=vtrace[du_idx,2][first_l0:last_l0] + vtracet=np.arange(0,len(vtracez))*dt_ns_l0[du_idx]+first + print("voltage",len(vtracez)) + + + atracex=atrace[du_idx,0][first_l1:last_l1]*109.86 #this number comes from doing 0.9V/8192, the maximum of ADC divided by 13bits + atracey=atrace[du_idx,1][first_l1:last_l1]*109.86 + atracez=atrace[du_idx,2][first_l1:last_l1]*109.86 + atracet=np.arange(0,len(atracez))*dt_ns_l1[du_idx]+first + print("adc",len(atracez)) + + + rtracex=rtrace[du_idx,0][first_l1:last_l1]*10 + rtracey=rtrace[du_idx,1][first_l1:last_l1]*10 + rtracez=rtrace[du_idx,2][first_l1:last_l1]*10 + rtracet=np.arange(0,len(rtracez))*dt_ns_l1[du_idx]+first + print("refield",len(rtracez)) fig,ax=plotfigure(time1=tracet,signal1=tracex,srate1=1/dt_ns_l0[du_idx],time2=vtracet,signal2=vtracex,srate2=1/dt_ns_l0[du_idx],time3=atracet,signal3=atracex,srate3=1/dt_ns_l1[du_idx],time4=rtracet,signal4=rtracex,srate4=1/dt_ns_l1[du_idx],label1="efield_l0 x10",label2="voltage_l0",label3="adc_l1 x110",label4="efield_l1 x10",Freqlimit=0.5) + #plt.savefig(args.directory+"IllustrateSimPipe"+str(du_idx)+".png") plt.show() diff --git a/sim2root/Common/sim2root.py b/sim2root/Common/sim2root.py index 3205c8b9..5033228c 100755 --- a/sim2root/Common/sim2root.py +++ b/sim2root/Common/sim2root.py @@ -34,6 +34,43 @@ print("#################################################") +def forcetimewindow(trace, t_bin_size, CurrentTPre,CurrentTPost,DesiredTPre=0,DesiredTPost=0): + ############################################################################################################################ + # adjust the trace lenght to force the requested tpre and tpost + ########################################################################################################################### + #print("before:",np.shape(efield),CurrentTPre,CurrentTPost, (CurrentTPre+CurrentTPost)/RawEfield.t_bin_size) + + if ForcedTPre!=0: + DeltaTimePre=ForcedTPre-CurrentTPre + DeltaBinsPre=int(np.round(DeltaTimePre/t_bin_size)) + else: + DeltaBinsPre=0 + + if ForcedTPost!=0: + DeltaTimePost=ForcedTPost-CurrentTPost + DeltaBinsPost=int(np.round(DeltaTimePost/t_bin_size)) + else: + DeltaBinsPost=0 + + if DeltaBinsPre<0 : + efield=efield[-DeltaBinsPre:] + DeltaBinsPre=0 + print("We have to remove "+str(-DeltaBinsPre)+" bins at the start of the trace") + if DeltaBinsPost<0 : + efield=efield[:DeltaBinsPost] + print("We have to remove "+str(-DeltaBinsPost)+" bins at the end of the trace") + + if DeltaBinsPost>0 or DeltaBinsPre>0: + npad = ((DeltaBinsPre, DeltaBinsPost), (0 , 0)) + trace=np.pad(trace, npad, 'constant') #TODO. Here I could pad using the "linear_ramp" mode, and pad wit a value that slowly aproached 0. + print("We have to add "+str(DeltaBinsPre)+" bins at the start of the trace") + print("We have to add "+str(DeltaBinsPost)+" bins at the end of the trace") + + #At this point, the traces are ensued to have a length between ForcedTpost and a Forced + + + + def convert_date(date_str): # Convert input string to a struct_time object date_struct = time.strptime(date_str, "%Y-%m-%d")