Algorithm

ps[T1_,T2_,T_,mn_]:=1-Exp[-T1 T mn]-Exp[-T2 T mn]+Exp[-mn (T1 T+T2 T-T1 T2 (T^2))];
pe[T1_,T2_,T_,mn_]:=(1-Exp[-T1 T mn]-Exp[-(1-T1) T mn]+Exp[-T mn]+1-Exp[-T2 T mn]-Exp[-(1-T2) T mn]+Exp[-T mn])/2
getPe[T1_,T2_,T_,qs_]:=Module[{i,step,mn},
i=0;mn=(10^(-8));step=0.4;
While[i<50&&step>0.000001,
If[(ps[T1,T2,T,mn]-qs)(ps[T1,T2,T,mn+step]-qs)<0,
step/=3,
mn+=step];
i++;

];
mn+=3step/2;
pe[T1,T2,T,mn]
]

See Fig. 1 in paper Quantum non-Gaussian coincidences

T1, T2 are transmission in the mode 1 and mode 2 respectively

T1=(Subscript[R, a,1])/((Subscript[R, a,1])+(Subscript[R, b,1]))

T2=(Subscript[R, a,2])/((Subscript[R, a,2])+(Subscript[R, b,2]))

(Subscript[R, a,1]) denotes click probability of (Subscript[SPAD, a,1]), (Subscript[R, b,1]) denotes click probability of (Subscript[SPAD, b,1]) etc.

T denotes transmission of the light before it reaches the detection part of the experiment, i. e. 1-T stands for the losses

Rs corresponds to success probability

T1=0.5;T2=0.5;T=1;Rs=1;
getPe[T1,T2,T,Rs]
Clear[T1];Clear[T2];Clear[T];Clear[Ts];Clear[Rs];

MapIndexed[Log10[{#1,getPe[1/2,1/2,0.5,#1]}]&,{0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.2}]

ListLinePlot[{{-2.`,-2.8755785947103227`},{-1.6989700043360187`,-2.336938649648952`},{-1.5228787452803376`,-2.037809968682972`},{-1.3979400086720375`,-1.8337124283299107`},{-1.3010299956639813`,-1.6802741057191322`},{-1.2218487496163564`,-1.5581009350755461`},{-1.154901959985743`,-1.4570474930218236`},{-1.0969100130080565`,-1.3711546403250343`},{-1.0457574905606752`,-1.296641393466932`},{-1.`,-1.2309671074816828`},{-0.958607314841775`,-1.1723406407892352`},{-0.9208187539523752`,-1.1194594740953459`},{-0.8860566476931632`,-1.0713444552703726`},{-0.8538719643217619`,-1.027245191021152`},{-0.8239087409443188`,-0.9865700673253592`},{-0.7958800173440752`,-0.9488494557570922`},{-0.769551078621726`,-0.913702213514349`},{-0.744727494896694`,-0.8808126313091147`},{-0.721246399047171`,-0.8499233928247931`},{-0.6989700043360187`,-0.820814376019068`}},AxesStyle->Directive[Large],ImageSize->600,AxesLabel->{"(Subscript[log, 10]) (Subscript[P, s])","(Subscript[log, 10]) ((Subscript[P, e1])+(Subscript[P, e2]))/2"},AxesOrigin->{-2,-3}]

(*this is a 4 detector scheme. the Rs is the probability of success of XX1-X1, T1(T2) is the BS transmission to XX1 (X1) and T are the losses before I reach the fiber BSs.*)

lrate=80000000;
blink=0.693598;
rabi=0.5;
drivingrate=lrate*blink*rabi

cc=219880/1200;
scc = Sqrt[219880]/1200.;
succ=cc/drivingrate
errcc=(348+285)/(1200*2);
serrcc=Sqrt[348^2+285^2]/(1200*2.);
rater=errcc/drivingrate

Analysis

data=Import["D:\\Checchinato\\NonGaussianity\\DATA_NONGAUSSIANITY\\rawdata\\Gaussianity_post_matlab_conv2_square11x11.txt","Data"];
data = Drop[data, 1];

Tsetup = 0.424;
Tdet = 0.5;
Tcoupling = 0.6919683881639405;
σTcoupling = 0.002182593342891458;
T1=Sqrt[0.5195725089060644*0.513736362870924];
T2=Sqrt[0.5788552646262161*0.5713013418201311];
T=Tsetup*Tdet*Tcoupling;
(*T=Sqrt[0.805161577205266*0.7795365625998585*Sqrt[0.5774423076712294*0.4384193141710676]*Sqrt[0.4615952846023682*0.5988890982383565]];*)
lrate=80000000;
blink=0.5660980066507492(*0.693598*);
σblink = 0.006212197431842366;
rabi=0.5;
drivingrate=lrate*blink*rabi*0.5;
intT = 1200;  
results = MapIndexed[{#[[1]],#[[2]],#[[3]]/(drivingrate*intT),(#[[4]]+#[[5]])/(intT*2*drivingrate),getPe[T1,T2,T,#[[3]]/(drivingrate*intT)],(#[[4]]+#[[5]])/(intT*2*drivingrate)/getPe[T1,T2,T,#[[3]]/(drivingrate*intT)] ,Abs[((#[[4]]+#[[5]])/(intT*2*drivingrate)-getPe[T1,T2,T,#[[3]]/(drivingrate*intT)])]/((#[[4]]+#[[5]])/(intT*2*drivingrate)*Sqrt[1/((#[[4]]+#[[5]])/(intT*2)*intT*2)+(σblink/blink)^2]),(#[[4]]+#[[5]] - (#[[6]]+#[[7]] ))/(intT*2*drivingrate)/getPe[T1,T2,T,#[[3]]/(drivingrate*intT)]  }&,data];
Clear[T1];Clear[T2];Clear[T];Clear[Ts];Clear[Rs];

TableForm[SortBy[results,#[[6]]&],TableHeadings->{None,{"tx(ps)","ty(ps)","success rate","error rate exp","error rate Theory","errExp/errTeo", "distance std", "Estimation no Background - Gaussianity_Jan.mat"}}]

Uncertainty Calculation

Uncertainty of experimental error/success rate : rate*Sqrt[(DeltaN/N)**2+(Delta_blink/blink)**2]

Uncertainty success: Sqrt[219880].;

Uncertainty experimental error: Sqrt[348^2+285^2]/2.;

Uncertainty (1 sigma) of  T1, T2, T, blinking

T1: 0.51665(26)

T2: 0.57508(35)

T: 0.424(16)

blinking: 0.6936(56)

Error estimation of  theoretical error using minimum and maximum values (3 sigma)

T1=0.51665 +σT1*3;
T2=0.57508 + σT2*3;
T=(Tsetup*Tdet*Tcoupling+ 3*σT);
Rs=succ - 3* succ*Sqrt[1/(cc*intT)+(σblink/blink)^2];
errmin = getPe[T1,T2,T,Rs];
Clear[T1];Clear[T2];Clear[T];Clear[Ts];Clear[Rs];

T1=0.51665 - σT1*3;
T2=0.57508 -σT2*3;
T=(Tsetup*Tdet*Tcoupling- 3*σT);
Rs=succ + 3* succ*Sqrt[1/(cc*intT)+(σblink/blink)^2];
errmax = getPe[T1,T2,T,Rs];
Clear[T1];Clear[T2];Clear[T];Clear[Ts];Clear[Rs];

T1=0.51665 ;
T2=0.57508 ;
T=0.5*0.424;
Rs=succ ;
errPe = getPe[T1,T2,T,Rs];
Clear[T1];Clear[T2];Clear[T];Clear[Ts];Clear[Rs];

succmax = 0.000017967516014016228`;
errormax =2.925730139554817`*^-8;
lrate=80000000;
(*blink=0.5721948156679171*)(*0.693598*)
rabi=0.5;
drivingrate=lrate*blink*rabi*0.5
σblink = 0.006212197431842366

cc = drivingrate *succmax;
(*cc=219880/1200;*)
succ=cc/drivingrate
σsucc = succ*Sqrt[1/(succ*drivingrate*1200)+(σblink/blink)^2];
errcc = drivingrate*errormax;
(*errcc=(348+285)/(1200*2);*)
rater=errcc/drivingrate
σrater = rater*Sqrt[1/(rater*drivingrate*1200*2)+(σblink/blink)^2];

σPe = ((errmax -errmin)/6)
σrater
rater/errPe*Sqrt[(σrater/rater)^2+(σPe/errPe)^2];
clear[rater,errPe,succ,errcc];

Tsetup = 0.424;
Tdet = 0.5;
Tcoupling = 0.6919683881639405;
σTcoupling = 0.002182593342891458;
T1=Sqrt[0.5195725089060644*0.513736362870924];
σT1 = 0.00025629615047562267;
T2=Sqrt[0.5788552646262161*0.5713013418201311];
σT2 = 0.00034834174650329576;
T=Tsetup*Tdet*Tcoupling;
σT = 0.0033658233299105156;
criterion = MapIndexed[{#1,getPe[T1,T2,T,#1]}&,10^Range[-1,1,0.01]*succmax];
criterionlowerbound = MapIndexed[{#1,getPe[T1,T2,T,#1]-(getPe[T1,T2,T,#1]-getPe[T1+3*σT1,T2+3*σT2,T+3*σT,#1*(1-3*Sqrt[1/(#1*drivingrate*1200)+(σblink/blink)^2])])/3}&,10^Range[-1,1,0.01]*succmax];
criterionupperbound = MapIndexed[{#1,getPe[T1,T2,T,#1]+(getPe[T1-3*σT1,T2-3*σT2,T-3*σT,#1*(1+3*Sqrt[1/(#1*drivingrate*1200)+(σblink/blink)^2])]-getPe[T1,T2,T,#1])/3}&,10^Range[-1,1,0.01]*succmax];
clear[T1,T2,T];

ListLinePlot[{criterionupperbound,criterion,criterionlowerbound, {Log10[{0.000017967516014016228`,2.925730139554817`*^-8}]}},AxesStyle->Directive[Large],ImageSize->600,AxesLabel->{"(Subscript[log, 10]) (Subscript[P, s])","(Subscript[log, 10]) ((Subscript[P, e1])+(Subscript[P, e2]))/2"},AxesOrigin->{-6,-10}]

Export["D:\\Checchinato\\NonGaussianity\\DATA_NONGAUSSIANITY\\Codes\\CodeOutput\\criterion.dat",criterion];
Export["D:\\Checchinato\\NonGaussianity\\DATA_NONGAUSSIANITY\\Codes\\CodeOutput\\criterionlowerbound.dat",criterionlowerbound];
Export["D:\\Checchinato\\NonGaussianity\\DATA_NONGAUSSIANITY\\Codes\\CodeOutput\\criterionupperbound.dat",criterionupperbound];

results = MapIndexed[{#[[1]]/(drivingrate*intT),(#[[2]]+#[[3]])/(intT*2*drivingrate),getPe[T1,T2,T,#[[1]]/(drivingrate*intT)],(#[[2]]+#[[3]])/(intT*2*drivingrate)/getPe[T1,T2,T,#[[1]]/(drivingrate*intT)] ,Abs[((#[[2]]+#[[3]])/(intT*2*drivingrate)-getPe[T1,T2,T,#[[1]]/(drivingrate*intT)])]/((#[[2]]+#[[3]])/(intT*2*drivingrate)*Sqrt[1/((#[[2]]+#[[3]])/(intT*2)*intT*2)+(σblink/blink)^2]),(#[[1]]/(drivingrate*1200))*Sqrt[1/((#[[1]]/(drivingrate*1200))*drivingrate*1200)+(σblink/blink)^2],(#[[2]]+#[[3]])/(intT*2*drivingrate)*Sqrt[1/((#[[2]]+#[[3]])/(intT*2*drivingrate)*drivingrate*1200*2)+(σblink/blink)^2]  }&,{{244113,365,430},{197776,308,348},{156550,240 , 275},{ 120180,180, 212}}];

TableForm[SortBy[results,#[[4]]&],TableHeadings->{None,{"success rate","error rate exp","error rate Theory","errExp/errTeo", "distance std", "σ success rate", "σ error rate exp"}}]

ATdata=Import["D:\\Checchinato\\NonGaussianity\\DATA_NONGAUSSIANITY\\rawdata\\FullTransmissionData.dat","Data"];
ATdata = Drop[ATdata, 1];
For [i=1,i<Length[ATdata]+1,i++,
AT = ATdata[[i,4]]/100.;
succmax = 0.000017967516014016228`;
errormax =2.925730139554817`*^-8;
Tsetup = 0.424;
Tdet = 0.5;
Tcoupling = 0.6919683881639405;
σTcoupling = 0.002182593342891458;
T1=Sqrt[0.5195725089060644*0.513736362870924];
σT1 = 0.00025629615047562267;
T2=Sqrt[0.5788552646262161*0.5713013418201311];
σT2 = 0.00034834174650329576;
T=Tsetup*Tdet*Tcoupling*AT;
σT = 0.0033658233299105156*AT;
criterion = MapIndexed[{#1,getPe[T1,T2,T,#1], getPe[T1,T2,T,#1]-(getPe[T1,T2,T,#1]-getPe[T1+3*σT1,T2+3*σT2,T+3*σT,#1])/3, getPe[T1,T2,T,#1]+(getPe[T1-3*σT1,T2-3*σT2,T-3*σT,#1]-getPe[T1,T2,T,#1])/3}&,10^Range[-2.4,1,0.01]*succmax];
Export[StringJoin[{"D:\\Checchinato\\NonGaussianity\\DATA_NONGAUSSIANITY\\Codes\\CodeOutput\\result",ToString[AT],".dat"}],Join[{{"success", "bound", "lowerbound", "upperbound"}},criterion]];
results = MapIndexed[{#[[4]],#[[1]]/(drivingrate*intT),(#[[2]]+#[[3]])/(intT*2*drivingrate),getPe[T1,T2,T,#[[1]]/(drivingrate*intT)],(#[[2]]+#[[3]])/(intT*2*drivingrate)/getPe[T1,T2,T,#[[1]]/(drivingrate*intT)] ,Abs[((#[[2]]+#[[3]])/(intT*2*drivingrate)-getPe[T1,T2,T,#[[1]]/(drivingrate*intT)])]/((#[[2]]+#[[3]])/(intT*2*drivingrate)*Sqrt[1/((#[[2]]+#[[3]])/(intT*2)*intT*2)+(σblink/blink)^2]),(#[[1]]/(drivingrate*1200))*Sqrt[1/((#[[1]]/(drivingrate*1200))*drivingrate*1200)+(σblink/blink)^2],(#[[2]]+#[[3]])/(intT*2*drivingrate)*Sqrt[1/((#[[2]]+#[[3]])/(intT*2*drivingrate)*drivingrate*1200*2)+(σblink/blink)^2]  }&,ATdata];
Export[StringJoin[{"D:\\Checchinato\\NonGaussianity\\DATA_NONGAUSSIANITY\\Codes\\CodeOutput\\result",ToString[AT],".dat"}],Join[{{"Simulated Transmission","success rate","error rate exp","error rate Theory","errExp/errTeo", "distance std", "σ success rate", "σ error rate exp"}},results]];
clear[T1,T2,T];
]

ATdata=Import["D:\\Checchinato\\NonGaussianity\\DATA_NONGAUSSIANITY\\rawdata\\FullTransmissionData.dat","Data"];

ATdata