debugged normalization in efield2voltage and efield2efield

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) 

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")

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()     

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")