<< NumericalMath`ListIntegrate`
<< Statistics`DataSmoothing`
<< NumericalMath`PolynomialFit`
<< Statistics`NormalDistribution`
<< Statistics`DiscreteDistributions`
hide = DisplayFunction -> Identity
SetDirectory["/afs/cern.ch/user/a/afabich/public/MERIT"]
<< "prolog.m"

\(\(\( (*\ 
      PS\ beam\ information\ from\ CPS\ *) \)\(
\)\(
\)\(\(PSbct = 
        ReadList["/afs/cern.ch/user/a/afabich/public/MERIT/BeamLog_CompactForm.txt", String];\)

    PSbct[\([1]\)]

    \(PSbct = 
        Transpose[Map[StringSplit, Drop[PSbct, 1]]];\)

    \(PSbct = 
        ToExpression[{PSbct[\([1]\)], PSbct[\([\(-3\)]\)], 
            PSbct[\([\(-2\)]\)], PSbct[\([\(-1\)]\)]}];\)

    \)\)\)

"Every shot with index smaller than 4014 has a two times lower \
attenuation, the signal is twice higher for the same intensity."
corr[x_] := 
      Join[Take[x, {1, 8}], 
        Take[x, {9, 24}]/If[1042 < x[[1]] < 4014, 2, 1], 
        Take[x, {25, 41}]]
repH8 = 
        ReadList["/afs/cern.ch/user/a/afabich/public/MERIT/transformer8.dat\
", Table[Number, {41}]];
repH8 = Map[corr, repH8];
repH16 = 
        Select[ReadList["/afs/cern.ch/user/a/afabich/public/MERIT/\
transformerWithSaturation.dat", Table[Number, {41}]], 
          IntegerPart[#[[1]]/1000] ≠ 
              20 &];
repH16 = Map[corr, repH16];
Export["/afs/cern.ch/user/a/afabich/public/MERIT/transformerAll.dat\
", report = Sort[Join[repH8, repH16]]]
report = Select[report, #[[-1]] == 0 &];

TotBctH16 = 
        Transpose[{Transpose[repH16][[1]], 
            Map[Total, 
              Transpose[
                Take[Transpose[repH16], {9, 
                    24}]]]}];
TotBctH16 = 
        Transpose[{Transpose[report][[1]], 
            Map[Total, 
              Transpose[Take[Transpose[report], {9, 24}]]]}];

sel[x_] := 
      Flatten[Join[
          Select[Transpose[
              PSbct], #[[1]] == 
                x[[1]] &], {x[[2]]}]]

    result = 
      Select[PSbctSel = Map[sel, TotBctH16], Length[#] == 5 &];   (* 
      remove runs with no data in PS bct *) 
result = Select[result, #[[5]]/#[[3]] > 100 &];   (* 
      remove runs with strange MERIT BCT data *) 
result = 
      Select[result, #[[1]] > 1999 &];  (* 
      remove runs earlier then 2000 *)

indexFull = 
        Transpose[{Table[
              m, {m, Length[Transpose[
                      result][[1]]]}]
, Transpose[
                result][[1]]}];
indexRed = {{indexFull[[1, 1]], 
            IntegerPart[indexFull[[1, 2]]/1000]}};
Do[
If[
        IntegerPart[indexFull[[m, 2]]/1000] - 
            indexRed[[-1, 2]] ≥ 1, 
        AppendTo[
          indexRed, {indexFull[[m, 1]], 
            IntegerPart[
              indexFull[[m, 2]]/1000]}]]
, {m, 2, 
        Length[indexFull]}]
indexRed

buf = Transpose[result][[5]]/Transpose[
                result][[2]]/
          Mean[Transpose[result][[5]]/Transpose[
                  result][[3]]];
buf2 = Transpose[result][[5]]/Transpose[
                result][[3]]/
          Mean[Transpose[result][[5]]/Transpose[
                  result][[3]]];
sel[x_] := 
      If[Abs[x[[1]] - 1] < Abs[x[[2]] - 1], x[[1]], 
        x[[1]]]
bufX = 
        Map[sel, Transpose[{buf, buf2}]]/0.904;

\(Export["/afs/cern.ch/user/a/afabich/public/MERIT/PicDisplayAllRuns.eps", 
Show[
        GraphicsArray[{
Histogram[bufX, 
              HistogramCategories -> 
                Table[m, {m, Floor[Min[bufX] 10]/10, 

                    Ceiling[Max[bufX] 10]/10, 0.005}], Frame -> True, 
              FrameLabel -> {"ratio MCT/PSCT", "entries"}, 
              hide], 
MultipleListPlot[
bufX, 
\(\(Transpose[result]\)[\([2]\)]/
                  Mean[\(Transpose[result]\)[\([2]\)]]\)/
                2, 
\(\(Transpose[result]\)[\([3]\)]/
                  Mean[\(Transpose[result]\)[\([3]\)]]\)/
                2, 
PlotRange -> {All, {0.0, 2}}, 
              FrameTicks -> {indexRed, True, False, False}, hide, 
              FrameLabel -> {"Run index (day)", "ratio MCT/PSCT"}, PlotJoined -> {False, True, True}, 
              SymbolShape -> {PlotSymbol[Box], None, None}, 
              PlotStyle -> {RGBColor[0, 0, 0], RGBColor[1, 0, 0], 
                  RGBColor[0, 1, 0], RGBColor[0, 1, 1], 
                  RGBColor[0, 0, 1]}]
}], 
        ImageSize -> 1100]
]\)

rzoom = {11000, 12999};
zoom = 
        Select[Transpose[
            Join[{Transpose[indexFull][[1]]}, {bufX}, {Transpose[
                    indexFull][[2]]}]], 
          rzoom[[1]] ≤ #[[3]] ≤ 
              rzoom[[2]] &];
prange = 
        Position[Transpose[
              indexFull][[2]], _?(rzoom[[1]] ≤ # ≤ rzoom[[2]] &)];
prange = 
        Join[First[prange], Last[prange]];

    
\(Export["/afs/cern.ch/user/a/afabich/public/MERIT/PicDisplayAllRunsZoom.eps", 
DisplayTogether[
ListPlot[
\(Transpose[
              zoom]\)[\([2]\)], 
Axes -> False, 
          Frame -> True, 
          FrameTicks -> {False, True, 
              False, {{0, 0}, {0.5, 3}, {1, 6}, {1.5, 9}}}, 
          PlotJoined -> False, 
          FrameLabel -> {"run index", "ratio MCT/PSCT (average normalized)", "", \*""total proton intensity [(10^12)]""}
], 
MultipleListPlot[
Take[\(Transpose[
                  result]\)[\([2]\)]/600, 
            prange], 
Take[\(Transpose[result]\)[\([3]\)]/
              600, prange]
, PlotJoined -> False, 
          SymbolShape -> PlotSymbol[Triangle], 
          SymbolStyle -> {RGBColor[1, 0, 0], 
              RGBColor[0, 1, 0]}], 
Show[
          Graphics[
Join[
              Table[Text[
                  zoom[\([m, 3]\)], {m, 
                    zoom[\([m, 2]\)]}, {0, \(-1\)}, {0, \(-1\)}], {m, 1, 
                  Length[zoom]}], {Text["intensity 3", {0, 0}, {\(-1\), 
                    0}, TextStyle -> {FontColor -> 
                        RGBColor[0, 1, 0]}]}, {Text["intensity 2", {0, 
                    0.1}, {\(-1\), 0}, 
                  TextStyle -> {FontColor -> 
                        RGBColor[1, 0, 0]}]}]
], 
          hide]
, ImageSize -> 700], 
      ImageRotated -> True]\)