## Analysis of Christensen et al. experiment using Graft's data-files
##
## Christensen, B.G., McCusker, K.T., Altepeter, J.B., Calkins, B., 
## Lim, C.C.W., Gisin, N., and Kwiat, P.G. (2013),
## "Detection-Loophole-Free Test of Quantum Nonlocality, and Applications", 
## Phys. Rev. Lett. 111, 130406.

## Donald A. Graft (2014)
## "Analysis of the Christensen et al. Clauser-Horne 
##      (CH)-Inequality-Based Test of Local Realism",
## http://arxiv.org/abs/1409.5158
## https://pubpeer.com/publications/E0F8384FC19A6034E86D516D03BB38



## Step 3: compute normalized Clauser-Horne B and B'




load("results.RData")


## Total counts, cf. Table 1 of Christensen et al. (2013)

result <- apply(counts, c(1, 2), sum)
t(result[c("SA", "C", "SB", "N"), ])
##        count
## setting     SA     C     SB        N
##      11  46070 29173  46202 27153018
##      12  48077 34146 146243 28352345
##      21 150837 34473  47448 27827311
##      22 150713  1869 144371 27926988
## Now a separate analysis of each run (sub-experiment)


for (exp in 1:20){
    result <- counts[ , , exp, drop = TRUE]

    N11 <- result["N", "11"]
    C11 <- result["C", "11"]
    S11A <- result["SA", "11"]
    S11B <- result["SB", "11"]

    N12 <- result["N", "12"]
    C12 <- result["C", "12"]
    S12A <- result["SA", "12"]
    S12B <- result["SB", "12"]

    N21 <- result["N", "21"]   
    C21 <- result["C", "21"]
    S21A <- result["SA", "21"]
    S21B <- result["SB", "21"]

    N22 <- result["N", "22"]   
    C22 <- result["C", "22"]
    S22A <- result["SA", "22"]
    S22B <- result["SB", "22"]

    S1dotA <- S11A + S12A
    Sdot1B <- S11B + S21B
    
    N1dot <- N11 + N12
    Ndot1 <- N11 + N21

    B <- (C11 / N11 + C12 / N12 + C21 / N21) - 
        (C22 / N22 + S1dotA / N1dot +Sdot1B / Ndot1) 
    
    Bprime <- (C11 / N11 + C12 / N12 + C21 / N21 - C22 / N22) / 
                    (S1dotA / N1dot + Sdot1B / Ndot1)
                    

    probs11 <- c(C11, S11A - C11, S11B - C11) / N11
    coefs11 <- c(1 - N11 / N1dot - N11 / Ndot1, - N11 / N1dot, - N11 / Ndot1)
    covmat11 <- diag(probs11) - outer(probs11, probs11)
    var11 <- sum(covmat11 * outer(coefs11, coefs11)) / N11

    probs12 <- c(C12, S12A - C12) / N12
    coefs12 <- c(1 - N12 / N1dot, - N12 / N1dot)
    covmat12 <- diag(probs12) - outer(probs12, probs12)
    var12 <- sum(covmat12 * outer(coefs12, coefs12)) / N12


    probs21 <- c(C21, S21B - C21) / N21
    coefs21 <- c(1 - N21 / Ndot1, - N21 / Ndot1)
    covmat21 <- diag(probs21) - outer(probs21, probs21)
    var21 <- sum(covmat21 * outer(coefs21, coefs21)) / N21

    probs22 <- C22 / N22
    var22 <- probs22 * (1 - probs22) / N22

    seB <- sqrt(var11 + var12 + var21 + var22)

    B2 <- sum(probs11 * coefs11) + sum(probs12 * coefs12) + sum(probs21 * coefs21) - probs22
    
    cat("\n", exp, B, seB, B/seB, Bprime)               
}
## 
##  1 7.693421e-05 5.688486e-05 1.352455 1.023532
##  2 0.0002450607 7.595572e-05 3.226363 1.074909
##  3 -1.172113e-06 6.144997e-05 -0.01907426 0.999638
##  4 9.56839e-05 1.957665e-05 4.887654 1.02834
##  5 -3.793885e-05 3.621009e-05 -1.047742 0.9892583
##  6 4.834904e-05 2.870192e-05 1.684523 1.013973
##  7 7.21832e-05 2.262471e-05 3.190458 1.021118
##  8 4.139045e-05 3.151019e-05 1.313557 1.01225
##  9 3.679404e-05 2.842475e-05 1.294437 1.010616
##  10 9.327091e-05 3.075131e-05 3.033071 1.027717
##  11 2.531185e-05 3.176529e-05 0.7968398 1.007528
##  12 2.118077e-05 2.982078e-05 0.7102688 1.006305
##  13 7.148607e-05 3.169182e-05 2.255663 1.021411
##  14 4.443709e-05 3.466292e-05 1.281978 1.013132
##  15 2.476795e-05 2.573426e-05 0.9624504 1.007367
##  16 4.369174e-05 2.252409e-05 1.939778 1.012766
##  17 6.098904e-05 2.722244e-05 2.240395 1.017868
##  18 6.30766e-05 3.056864e-05 2.063442 1.018499
##  19 4.762152e-05 2.982972e-05 1.596446 1.013618
##  20 3.563453e-05 3.190108e-05 1.117032 1.01064
cat("\nDataSet   B        se(B)       B/se(B)       B'")
## 
## DataSet   B        se(B)       B/se(B)       B'
## Experiments 3 and 5 are anomalous.
## Here are the counts for these two experiments.

## Z is the number of windows with no events at all.
## N is the total number of windows.
## C is the number of windows with coincidences (detections by Alice and by Bob).
## SA and SB are the numbers of windows with detections by Alice and Bob respectively
## (so C is included in both SA and SB). 

## So we should have N = Z + SA + SB - C


counts[ , , 3]
##      setting
## count     11     12     21     22
##    Z  199578 248633 372935 420729
##    N  200028 250038 375050 425011
##    C     192    285    452     29
##    SA    329    409   1961   2245
##    SB    313   1281    606   2066
counts[ , , 5]
##      setting
## count      11     12      21      22
##    Z  1447101 919834 1292416 1311232
##    N  1450827 925066 1300120 1325132
##    C     1566   1166    1581      82
##    SA    2480   1626    7135    7087
##    SB    2812   4772    2150    6895