## 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 2: just get all coincidence counts and singles counts and empty window counts.
## Save all counts of all experiments in one binary file


summary <- read.csv("summary.csv", stringsAsFactors = FALSE)
Names <- summary$Name

lengths <- summary$length



windowDiamMS <- 2.4 * 10^{-6}
delayA <- (2.65 * 10^{-6} - windowDiamMS/2) / (156.25 * 10^{-12})
delayB <- (2.525 * 10^{-6} - windowDiamMS/2) / (156.25 * 10^{-12})

delayA <- as.integer(round(delayA))
delayB <- as.integer(round(delayB))

Window <- windowDiamMS / (156.25 * 10^{-12})
Window <- as.integer(round(Window))

alice <- function(x) any(x == 1, na.rm = TRUE)
bob <- function(x) any(x == 2, na.rm = TRUE)

Nrep <- 20

results <- list()

for (Rep in 1:Nrep) {
    Name <- Names[Rep]  
    File <- paste(substr(Name, 1, 6), ".RData", sep ="")
    cat("\n", "\n", File)
    load(File)

    N <- length(Times)
    

    while (Events[1] != 0) {
        Times <- Times[-1]
        Events <- Events[-1]
        Settings <- Settings[-1]
        N <- N - 1
    }
    
    while (Events[N] == 0) {
        Times <- Times[-N]
        Events <- Events[-N]
        Settings <- Settings[-N]
        N <- N - 1
    }    
        

    i <- 0
    Count <- 0
    maxEvents <- 0
    nBunch <- 1
    t <- 0
    while (i < N) {
      i <- i + 1
      if (Events[i] == 0) {
        maxEvents <- max(nBunch, maxEvents)
        nBunch <- 1
        Count <- Count + 1
        t <- Diffs[i]
      } else {
        t <- t + Diffs[i]
        if (t < Window) nBunch <- nBunch + 1
      }
    }          
    maxEvents <- max(nBunch, maxEvents)
    
    cat("\n", Count, "x", maxEvents)

    NewEvents <- matrix(NA, Count, maxEvents)
    NewSettings <- integer(Count)


    i <- 0
    j <- 0
    k <- 0
    t <- 0
    while (i < N) {
      i <- i + 1
      if (Events[i] == 0) {
        j <- j + 1
        k <- 0
        t <- 0
        NewSettings[j] <- Settings[i]
      } else {
        t <- t + Diffs[i - 1]
        if (t < Window) {
          k <- k + 1
          NewEvents[j, k] <- Events[i]
        }
      }
    }          

    Alice <- apply(NewEvents, 1, alice)
    Bob <- apply(NewEvents, 1, bob)

    N11 <- nWindows[1]
    C11 <- sum((Alice * Bob)[NewSettings == 11])
    S11A <- sum(Alice[NewSettings == 11])
    S11B <- sum(Bob[NewSettings == 11])
    Z11 <- N11 - S11A - S11B + C11
    cat("\n", "11:", C11, S11A, S11B, Z11)

    N12 <- nWindows[2]
    C12 <- sum((Alice * Bob)[NewSettings == 12])
    S12A <- sum(Alice[NewSettings == 12])
    S12B <- sum(Bob[NewSettings == 12])
    Z12 <- N12 - S12A - S12B + C12
    cat("\n", "12:", C12, S12A, S12B, Z12)

    N21 <- nWindows[3]
    C21 <- sum((Alice * Bob)[NewSettings == 21])
    S21A <- sum(Alice[NewSettings == 21])
    S21B <- sum(Bob[NewSettings == 21])
    Z21 <- N21 - S21A - S21B + C21
    cat("\n", "21:", C21, S21A, S21B, Z21)

    N22 <- nWindows[4]
    C22 <- sum((Alice * Bob)[NewSettings == 22])
    S22A <- sum(Alice[NewSettings == 22])
    S22B <- sum(Bob[NewSettings == 22])
    Z22 <- N22 - S22A - S22B + C22
    cat("\n", "22:", C22, S22A, S22B, Z22)  
   
    results[[Rep]] <- list(list(Z11, N11, C11, S11A, S11B), 
                           list(Z12, N12, C12, S12A, S12B),
                           list(Z21, N21, C21, S21A, S21B),
                           list(Z22, N22, C22, S22A, S22B)
    )
                            
}
## 
##  
##  data01.RData
##  8019 x 5
##  11: 266 410 392 249492
##  12: 477 673 2023 397807
##  21: 385 1676 530 323216
##  22: 26 1422 1456 272172
##  
##  data02.RData
##  8184 x 5
##  11: 317 490 489 324377
##  12: 516 709 2069 397788
##  21: 230 965 320 173949
##  22: 22 1824 1747 346489
##  
##  data03.RData
##  9016 x 5
##  11: 192 329 313 199578
##  12: 285 409 1281 248633
##  21: 452 1961 606 372935
##  22: 29 2245 2066 420729
##  
##  data04.RData
##  83179 x 9
##  11: 3469 5485 5384 3217842
##  12: 3694 5151 15502 3058347
##  21: 4096 17769 5590 3255982
##  22: 162 15692 14911 2894765
##  
##  data05.RData
##  35169 x 13
##  11: 1566 2480 2812 1447101
##  12: 1166 1626 4772 919834
##  21: 1581 7135 2150 1292416
##  22: 82 7087 6895 1311232
##  
##  data06.RData
##  34915 x 6
##  11: 1105 1741 1763 997677
##  12: 1754 2495 7248 1417121
##  21: 1857 8084 2521 1491371
##  22: 70 5845 5521 1063734
##  
##  data07.RData
##  74024 x 10
##  11: 2566 4079 4000 2344719
##  12: 3179 4399 13534 2610502
##  21: 2781 12154 3800 2187001
##  22: 216 15612 14793 2794886
##  
##  data08.RData
##  34930 x 8
##  11: 1402 2258 2201 1322049
##  12: 1472 2098 6583 1242872
##  21: 1598 6900 2186 1267626
##  22: 79 6139 5999 1138027
##  
##  data09.RData
##  34398 x 13
##  11: 1222 1932 1874 1097498
##  12: 1794 2586 7778 1516554
##  21: 2063 8846 2838 1590521
##  22: 57 4301 4466 766353
##  
##  data10.RData
##  35629 x 6
##  11: 1256 1985 2004 1172383
##  12: 1803 2483 7585 1441762
##  21: 1323 5918 1780 1093717
##  22: 72 6875 6682 1261655
##  
##  data11.RData
##  36697 x 6
##  11: 1195 1908 1864 1122533
##  12: 1236 1722 5330 1019257
##  21: 1749 7849 2447 1441593
##  22: 121 7657 7279 1385257
##  
##  data12.RData
##  33334 x 8
##  11: 1412 2176 2226 1297100
##  12: 1759 2477 7678 1466691
##  21: 1530 6797 2150 1292719
##  22: 69 5102 4732 915299
##  
##  data13.RData
##  35344 x 8
##  11: 1322 2108 2139 1272209
##  12: 1519 2113 6777 1292769
##  21: 1450 6238 1983 1143332
##  22: 71 6957 6454 1261724
##  
##  data14.RData
##  35496 x 7
##  11: 1443 2215 2320 1347000
##  12: 1248 1768 5355 1044219
##  21: 1399 6116 1948 1118401
##  22: 107 7792 7535 1459879
##  
##  data15.RData
##  52026 x 16
##  11: 2048 3290 3198 1945734
##  12: 2119 3051 9218 1764992
##  21: 2474 10801 3400 2013501
##  22: 138 9425 9098 1731794
##  
##  data16.RData
##  70881 x 8
##  11: 2730 4270 4273 2469411
##  12: 3364 4692 14146 2734767
##  21: 2627 11767 3667 2162331
##  22: 176 13856 13382 2573138
##  
##  data17.RData
##  54410 x 8
##  11: 2045 3214 3239 1895732
##  12: 2261 3233 9961 1914247
##  21: 1996 8404 2723 1540974
##  22: 136 11502 10864 2102924
##  
##  data18.RData
##  36709 x 8
##  11: 1235 1910 1921 1122490
##  12: 1502 2107 6482 1243067
##  21: 1530 6742 2161 1242703
##  22: 69 7273 7029 1360880
##  
##  data19.RData
##  37118 x 6
##  11: 1155 1800 1797 1022656
##  12: 1516 2141 6459 1218006
##  21: 1874 8233 2611 1491137
##  22: 86 6844 6501 1236845
##  
##  data20.RData
##  36267 x 9
##  11: 1227 1990 1993 1222338
##  12: 1482 2144 6462 1242936
##  21: 1478 6482 2037 1168074
##  22: 81 7263 6961 1335991
counts <- array(0, dim=c(5, 4, 20), 
                dimnames = list(count = c("Z", "N", "C", "SA", "SB"), 
                                setting = c("11", "12", "21", "22"),
                                experiment = 1:20
                                )
                )             
                
for (exp in 1:20){ 
    for (set in 1:4){
        for (c in 1:5){
            counts[c, set, exp] <- results[[exp]][[set]][[c]]
        }
    }
}   


save(counts, file = "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