(* 17 May 2014: This Mathematica code can be used as a quick check of the results obtained by Hiroyuki's and Francesc's codes. *)

(* 17 Sept. 2013: instead of finding triple product and then summing, as in gammas-northsouth notebook, here we first sum and then take triple product. *)
(* 13 Sept. 2013: To calculate the "helicity statistic" for gamma ray data for the project with Hiroyuki and Francesc and Wenlei. *)

(* Import data files from which Q will be calculated. The gamma ray data is assumed to have the sources removed. Insert datafile addresses. *)
d10=Import["~/Desktop/data_set/LAT_10-20GeV.dat","Data"]; 
d20=Import["~/Desktop/data_set/LAT_20-30GeV.dat","Data"];
d30=Import["~/Desktop/data_set/LAT_30-40GeV.dat","Data"];
d40=Import["~/Desktop/data_set/LAT_40-50GeV.dat","Data"];
d50=Import["~/Desktop/data_set/LAT_50-60GeV.dat","Data"];
 (*length of each data set*)
l10=Length[d10];
l20=Length[d20];
l30=Length[d30];
l40=Length[d40];
l50=Length[d50];
(*ListPlot[Table[d50[[i]][[2]],{i,1,l50}]]*)
(* spherical theta angle data -- north pole *)
th10=Table[Pi/2-d10[[n,2]],{n,1,l10}];
th20=Table[Pi/2-d20[[n,2]],{n,1,l20}];
th30=Table[Pi/2-d30[[n,2]],{n,1,l30}];
th40=Table[Pi/2-d40[[n,2]],{n,1,l40}];
th50=Table[Pi/2-d50[[n,2]],{n,1,l50}];
(*azimuthal data*)
ph10=Table[d10[[n,1]],{n,1,l10}];
ph20=Table[d20[[n,1]],{n,1,l20}];
ph30=Table[d30[[n,1]],{n,1,l30}];
ph40=Table[d40[[n,1]],{n,1,l40}];
ph50=Table[d50[[n,1]],{n,1,l50}];
(*construct vectors on the sphere*)v10=Table[{Sin[th10[[n]]]*Cos[ph10[[n]]],Sin[th10[[n]]]*Sin[ph10[[n]]],Cos[th10[[n]]]},{n,1,l10}];
v20=Table[{Sin[th20[[n]]]*Cos[ph20[[n]]],Sin[th20[[n]]]*Sin[ph20[[n]]],Cos[th20[[n]]]},{n,1,l20}];
v30=Table[{Sin[th30[[n]]]*Cos[ph30[[n]]],Sin[th30[[n]]]*Sin[ph30[[n]]],Cos[th30[[n]]]},{n,1,l30}];
v40=Table[{Sin[th40[[n]]]*Cos[ph40[[n]]],Sin[th40[[n]]]*Sin[ph40[[n]]],Cos[th40[[n]]]},{n,1,l40}];
v50=Table[{Sin[th50[[n]]]*Cos[ph50[[n]]],Sin[th50[[n]]]*Sin[ph50[[n]]],Cos[th50[[n]]]},{n,1,l50}];
(*datadist=Table[v50[[i]][[3]],{i,1,l50}];
ListPlot[datadist]
Histogram[datadist,100]*)
(* We only want patches higher than a certain latitude. So exclude 50 GeV photons close to the equator. *)
exclude=80; (* angle to be excluded *)
edgeplus=Cos[Pi/2.-exclude*Pi/180.];
e5=Select[v50,Abs[#[[3]]]>edgeplus&];
le5=Length[e5](* number of 50 GeV patches *)

(* Initialize arrays for storing Q and its standard error. *)
qval={};
sval={};

(* this is for 10, 20, 50 GeV photons, for patch radius = "patch". *)
patch=10;(* patch radius in degrees *)
pc=Cos[patch*Pi/180.] ;
(* select 20 GeV photons inside patch, for every patch. *)
e25=Table[Select[v20,#.e5[[i]]>pc&],{i,1,le5}];
(* select 10 GeV photons inside patch, for every patch. *)
e15=Table[Select[v10,#.e5[[i]]>pc&],{i,1,le5}];
(* sum over 10 GeV photons inside patch, for all patches. *)
pe15=Table[Sum[ e15[[k]][[i]],{i,1,Length[e15[[k]]]}],{k,1,Length[e5]}];
(* sum over 20 GeV photons inside patch, for all patches. *)
pe25=Table[Sum[ e25[[k]][[j]],{j,1,Length[e25[[k]]]}],{k,1,Length[e5]}];
(* evaluate triple product in every patch. *)
hinpatches125=Table[pe15[[k]]\[Cross]pe25[[k]].e5[[k]]/(Length[e15[[k]]]*Length[e25[[k]]]),{k,1,Length[e5]}];

(* Q value for (10,20,50) GeV; spread; standard error *)
stat125=Mean[hinpatches125] (* Q statistics *)
spread125=StandardDeviation[hinpatches125];
sdofmean125=spread125/Sqrt[Length[e5]]
(*Histogram[hinpatches125];*)

(* this section for plots. Have to run above routine for various patch radii. *)

PrependTo[qval,stat125]
PrependTo[sval,sdofmean125] (* first value in tables is for patch=2 *)

ListPlot[{qval,sval}]

plotData=Partition[Riffle[qval,sval],2]

Needs["ErrorBarPlots`"]

tickSpecs={{{1,"2"},{2,"3"},{3,"4"},{4,"5"},{5,"6"},{6,"7"},{7,"8"},{8,"9"},{9,"10"},{10,"11"},{11,"12"}},Automatic};

ErrorListPlot[plotData,Joined->True,Frame->{Left},FrameLabel->{None,Style["Q",12,FontFamily->"Arial"]}]

(* this is for 30, 40, 50 GeV photons *)
e45=Table[Select[v40,#.e5[[i]]>pc&],{i,1,le5}];
e35=Table[Select[v30,#.e5[[i]]>pc&],{i,1,le5}];
pe35=Table[Sum[ e35[[k]][[i]],{i,1,Length[e35[[k]]]}],{k,1,Length[e5]}];
pe45=Table[Sum[ e45[[k]][[j]],{j,1,Length[e45[[k]]]}],{k,1,Length[e5]}];
(*hinpatches345=Table[pe35[[k]]\[Cross]pe45[[k]].e5[[k]]/(Length[e35[[k]]]*Length[e45[[k]]]),{k,1,Length[e5]}];*)
hinpatches345=Table[If[Length[e35[[k]]]*Length[e45[[k]]] > 0,pe35[[k]]\[Cross]pe45[[k]].e5[[k]]/(Length[e35[[k]]]*Length[e45[[k]]]),0],{k,1,Length[e5]}];
patch
exclude
le5
stat345=Mean[hinpatches345] (* statistics evaluated for real data *)
spread345=StandardDeviation[hinpatches345];
sdofmean345=spread345/Sqrt[Length[e5]]
(*Histogram[hinpatches345];*)
qstat345[patch]=stat345;
sstat345[patch]=sdofmean345;