@@ -762,10 +762,10 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
762762 ln = '#avg_time_res, d_avg_time_res, var_time, var_ref, Correction\n '
763763 file_results .write (ln )
764764
765- best_result [k ] = Bauble ()
766- best_result [k ].dT = [- 1 , - 1 ]
767- best_result [k ].var = [None , None ]
768- best_result [k ].Correction = None
765+ best_result [N_bar ] = Bauble ()
766+ best_result [N_bar ].dT = [- 1 , - 1 ]
767+ best_result [N_bar ].var = [None , None ]
768+ best_result [N_bar ].Correction = None
769769
770770 if confL ['idx_ref' ] != confR ['idx_ref' ]:
771771 print 'Inconsistency in reference chnnel'
@@ -776,6 +776,7 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
776776
777777 time_var_chref = v_ref + '[{}]' .format (confL ['idx_ref' ])
778778 time_var = '0.5*({v}[{L}]+{v}[{R}])' .format (v = v_time , L = kL , R = kR )
779+ spacelike_var = '({v}[{L}]-{v}[{R}])' .format (v = v_time , L = kL , R = kR )
779780
780781 out_dir = '{}/{}_{}' .format (headout_dir , v_time , v_ref )
781782 if not os .path .exists (out_dir ):
@@ -809,6 +810,51 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
809810 print 'Not enought stat ({})' .format (chain .GetEntries (selection ))
810811 continue
811812
813+ '''=========================== Space like variable ==========================='''
814+ if 'SpaceLike' in configurations .plots :
815+ aux_x = tree2array (chain , confL ['x_rot' ], selection ).flatten ()
816+ aux_x -= np .min (aux_x )
817+ aux_spacelike = tree2array (chain , spacelike_var , selection ).flatten ()
818+
819+ inputs = np .column_stack ((np .ones_like (aux_x ), aux_x ))
820+ coeff , r , rank , s = np .linalg .lstsq (inputs , aux_spacelike , rcond = None )
821+ b = [None , None , None , None , np .percentile (aux_spacelike , 1 ), np .percentile (aux_spacelike , 99 )]
822+ h_SLvsX = create_TH2D (np .column_stack ((aux_x , aux_spacelike )),
823+ 'h_SLvsX' , 'Bar ' + str (N_bar ),
824+ binning = b ,
825+ axis_title = ['x [mm]' , spacelike_var + ' [ns]' ]
826+ )
827+ h_SLvsX .SetStats (0 )
828+ canvas ['dt_vs_amp' ][k ] = rt .TCanvas ('h_SLvsX' + str (N_bar ), 'h_SLvsX' + str (N_bar ), 800 , 600 )
829+ h_SLvsX .DrawCopy ('colz' )
830+ prof = h_SLvsX .ProfileX ('prof_amp' )
831+ prof .SetLineColor (6 )
832+ prof .SetLineWidth (2 )
833+ prof .DrawCopy ('SAMEE1' )
834+
835+ f = rt .TF1 ('SLvsX_fit' + str (k ),'[0]+[1]*x' , np .min (aux_x ), np .max (aux_x ))
836+ for j ,a in enumerate (coeff ):
837+ f .SetParameter (j , a )
838+ f .SetLineColor (6 )
839+ f .SetLineStyle (9 )
840+ f .DrawCopy ('SAMEL' )
841+
842+ l = rt .TLatex ()
843+ l .SetTextSize (0.04 )
844+ v = 2. / coeff [1 ]
845+ length = - coeff [0 ]* v
846+ v2 = - 50. / coeff [0 ]
847+ v3 = np .sum (np .square (2 * aux_x - 50. ))/ np .sum ((2 * aux_x - 50. )* aux_spacelike )
848+ print 'v3 = {:.0f} mm/ns' .format (v3 )
849+ msg = 'v = {:.0f} ({:.0f}) mm/ns - l = {:.1f} mm' .format (v , v2 , length )
850+ print msg
851+ print 'Refractive index = {:.2f}' .format (299.0 / v3 )
852+ l .DrawLatexNDC (0.15 , 0.89 , msg )
853+
854+ canvas ['dt_vs_amp' ][k ].Update ()
855+ canvas ['dt_vs_amp' ][k ].SaveAs (out_dir + '/SLvsX_{}.png' .format (N_bar ))
856+
857+
812858 '''=========================== Raw time resolution ==========================='''
813859 if 'TimeResRaw' in configurations .plots :
814860 delta_t = tree2array (chain , var_dT , selection ).flatten ()
@@ -848,16 +894,17 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
848894 print 'Please specify a shape for both ends in order to print TimeResRaw2D'
849895 continue
850896
851- b = [var_dT , conf ['x_rot' ], conf ['y_rot' ]]
897+ b = [var_dT , confL ['x_rot' ], confL ['y_rot' ]]
852898 if 'TimeRes2DAmp' in configurations .plots :
853899 b += ['amp[{}]' .format (confL ['amp_ch' ]), 'amp[{}]' .format (confR ['amp_ch' ])]
900+ b += [spacelike_var ]
854901
855902 data_raw = tree2array (chain , b , selection ).view (np .recarray )
856903 data = np .zeros ((data_raw .shape [0 ],3 ))
857904 if 'TimeRes2DAmp' in configurations .plots :
858- data = np .zeros ((data_raw .shape [0 ],5 ))
905+ data = np .zeros ((data_raw .shape [0 ],6 ))
859906 for iev , d in enumerate (data_raw ):
860- data [iev ] = [d [0 ][0 ], d [1 ], d [2 ], d [3 ], d [4 ]]
907+ data [iev ] = [d [0 ][0 ], d [1 ], d [2 ], d [3 ], d [4 ], d [ 5 ] ]
861908 else :
862909 for iev , d in enumerate (data_raw ):
863910 data [iev ] = [d [0 ][0 ], d [1 ], d [2 ]]
@@ -928,12 +975,16 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
928975
929976 if 'TimeRes2DAmp' in configurations .plots :
930977 dt_sel = np .logical_and (dt > median - 2 * width , dt < median + 2 * width )
931- amp = aux_d [:,3 ]/ aux_d [:,4 ]
932- inputs = np .column_stack ((np .ones_like (amp ), amp , amp ** 2 ))
978+ ampL = aux_d [:,3 ]
979+ ampR = aux_d [:,4 ]
980+ x_bar = aux_d [:,5 ]
981+ amp_ratio = ampL / ampR
982+ # inputs = np.column_stack((np.ones_like(amp_ratio), amp_ratio, amp_ratio**2))
983+ inputs = np .column_stack ((np .ones_like (ampL ), ampL , ampR , np .square (ampL ), ampR * ampL , np .square (ampR )))
933984 coeff , r , rank , s = np .linalg .lstsq (inputs [dt_sel ], dt [dt_sel ], rcond = None )
934985 b = [None , None , None , None , median - 2 * width , median + 2 * width ]
935- h_TvsA = create_TH2D (np .column_stack ((amp , dt )),
936- 'h_TvsA' , 'Channel ' + str (k ),
986+ h_TvsA = create_TH2D (np .column_stack ((amp_ratio , dt )),
987+ 'h_TvsA' , 'Bar ' + str (N_bar ),
937988 binning = b ,
938989 axis_title = ['Amp{} (L)/ Amp{} (R)' .format (kL , kR ), var_dT + ' [ns]' ]
939990 )
@@ -944,12 +995,12 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
944995 prof .SetLineWidth (2 )
945996 prof .DrawCopy ('SAMEE1' )
946997
947- f = rt .TF1 ('amp_fit' + str (k ),'[0]+[1]*x+[2]*x^2' , np .min (amp ), np .max (amp ))
948- for j ,a in enumerate (coeff ):
949- f .SetParameter (j , a )
950- f .SetLineColor (6 )
951- f .SetLineStyle (9 )
952- f .DrawCopy ('SAMEL' )
998+ # f = rt.TF1('amp_fit'+str(k),'[0]+[1]*x+[2]*x^2', np.min(amp), np.max(amp))
999+ # for j,a in enumerate(coeff):
1000+ # f.SetParameter(j, a)
1001+ # f.SetLineColor(6)
1002+ # f.SetLineStyle(9)
1003+ # f.DrawCopy('SAMEL')
9531004 canvas ['dt_vs_amp' ][k ].Update ()
9541005 canvas ['dt_vs_amp' ][k ].SaveAs (out_dir + '/TimeRes2DAmp_bar{:02d}' .format (N_bar )+ '/TvsAmp_{}.png' .format (ib ))
9551006
@@ -980,10 +1031,10 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
9801031 results .append ([avg_time_res , d_avg_time_res , v_time , v_ref , False ])
9811032 ln = '{:.2f} {:.2f} {} {} {}\n ' .format (avg_time_res , d_avg_time_res , v_time , v_ref , False )
9821033 file_results .write (ln )
983- if best_result [k ].dT [0 ] < 0 or best_result [k ].dT [0 ] > avg_time_res :
984- best_result [k ].dT = [avg_time_res , d_avg_time_res ]
985- best_result [k ].var = [v_time , v_ref ]
986- best_result [k ].AmpCorr = False
1034+ if best_result [N_bar ].dT [0 ] < 0 or best_result [N_bar ].dT [0 ] > avg_time_res :
1035+ best_result [N_bar ].dT = [avg_time_res , d_avg_time_res ]
1036+ best_result [N_bar ].var = [v_time , v_ref ]
1037+ best_result [N_bar ].AmpCorr = False
9871038
9881039 msg = 'Avg raw resolution: {:.1f} +/- {:.1f} ps' .format (avg_time_res , d_avg_time_res )
9891040 print msg
@@ -1019,33 +1070,22 @@ def fill_TimeResHisto(k, dt, h2D, out_list, tagin, title_tag, canvas, save_canva
10191070 results .append ([avg_time_res , d_avg_time_res , v_time , v_ref , True ])
10201071 ln = '{:.2f} {:.2f} {} {} {}\n ' .format (avg_time_res , d_avg_time_res , v_time , v_ref , True )
10211072 file_results .write (ln )
1022- if best_result [k ].dT [0 ] < 0 or best_result [k ].dT [0 ] > avg_time_res :
1023- best_result [k ].dT = [avg_time_res , d_avg_time_res ]
1024- best_result [k ].var = [v_time , v_ref ]
1025- best_result [k ].AmpCorr = True
1073+ if best_result [N_bar ].dT [0 ] < 0 or best_result [N_bar ].dT [0 ] > avg_time_res :
1074+ best_result [N_bar ].dT = [avg_time_res , d_avg_time_res ]
1075+ best_result [N_bar ].var = [v_time , v_ref ]
1076+ best_result [N_bar ].AmpCorr = True
10261077 msg = 'Avg resolution after amp correction: {:.1f} +/- {:.1f} ps' .format (avg_time_res , d_avg_time_res )
10271078 print msg
10281079 l = rt .TLatex ()
10291080 l .SetTextSize (0.04 );
10301081 l .DrawLatexNDC (0.15 , 0.89 , msg .replace ('+/-' , '#pm' ))
1031- if 'SiPM_sel' in conf .keys ():
1032- x_start , x_stop , y_start , y_stop = conf ['SiPM_sel' ]['limits' ]
1033- elif 'amp_channel' in conf .keys ():
1034- x_start , x_stop , y_start , y_stop = configurations .channel [conf ['amp_channel' ]]['SiPM_sel' ]['limits' ]
1035- line .SetLineStyle (9 )
1036- line .SetLineColor (46 )
1037- line .SetLineWidth (2 )
1038- line .DrawLine (x_start , y_start , x_stop , y_start )
1039- line .DrawLine (x_start , y_stop , x_stop , y_stop )
1040- line .DrawLine (x_start , y_start , x_start , y_stop )
1041- line .DrawLine (x_stop , y_start , x_stop , y_stop )
10421082 canvas ['c_' + tag ].Update ()
10431083 canvas ['c_' + tag ].SaveAs (out_dir + '/TimeRes2DAmp_bar{:02d}.png' .format (N_bar ))
10441084
10451085 file_results .close ()
10461086
10471087 print '\n \n ======================= Summary =============================='
1048- table = PrettyTable (['Ch ' , 'Best Resolution [ps]' , 'Var ref' , 'Var timr' , 'Amp corrected' ])
1088+ table = PrettyTable (['Bar ' , 'Best Resolution [ps]' , 'Var ref' , 'Var timr' , 'Amp corrected' ])
10491089 for k , res in best_result .iteritems ():
10501090 row = [str (k ), '{:.2f} +/- {:.2f}' .format (res .dT [0 ],res .dT [1 ])]
10511091 row += [res .var [1 ], res .var [0 ], 'Yes' if res .AmpCorr else 'No' ]
0 commit comments