#include <stdlib.h>
#include <fstream.h>
#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#include "Riostream.h"
#include "TH1F.h"
#include "TFile.h"
#include "TTree.h"
class Angle
{public:	Angle();
			double calc_angle(double, double, double, double);
};
Angle::Angle() { };
Angle::calc_angle(double ra1,double dec1,double ra2,double dec2)
{
double alp1,bet1,gam1,alp2,bet2,gam2,cost;
double rfactor=180./3.1415928;
alp1=cos(ra1/rfactor)*cos(dec1/rfactor);
bet1=sin(ra1/rfactor)*cos(dec1/rfactor);
gam1=sin(dec1/rfactor);
alp2=cos(ra2/rfactor)*cos(dec2/rfactor);
bet2=sin(ra2/rfactor)*cos(dec2/rfactor);
gam2=sin(dec2/rfactor);
cost=(alp1*alp2+bet1*bet2+gam1*gam2);
if(cost>1.0)cost=1.0;
cost=rfactor*acos(cost);
return cost;
};

void test_jdsummary(void)
{
gROOT->Reset();

short levent,leventf,levents;
float ra1,dec1,ra2,dec2,,delx;
UChar_t cutf,cuts,ntrackf,ntracks;
int runevent;
double jd,jdtemp,ramoon,decmoon,rasun,decsun,ra_m[240000],decl_m[240000];
double ra_s[10000],decl_s[10000],rasun,decsun,temp,moonangle,sunangle;
int nevent,i,j,index;
Angle a;

c1 = new TCanvas("c1","The Test Canvas",200,10,700,500);
TH1F *h1 = new TH1F(" moon angle ", " moon angle ",50,0.,5.);
TFile f("jdsummary.root");
TTree *newtree = (TTree*)f.Get("data");
data->SetBranchAddress("ra1",&ra1);
data->SetBranchAddress("dec1",&dec1);
data->SetBranchAddress("ra2",&ra2);
data->SetBranchAddress("dec2",&dec2);
data->SetBranchAddress("jd",&jd);
data->SetBranchAddress("cutf",&cutf);
data->SetBranchAddress("cuts",&cuts);
data->SetBranchAddress("levent",&levent);
data->SetBranchAddress("leventf",&leventf);
data->SetBranchAddress("levents",&levents);
data->SetBranchAddress("ntrackf",&ntrackf);
data->SetBranchAddress("ntracks",&ntracks);
data->SetBranchAddress("runevent",&runevent);


ifstream moondatafile("moon_position.dat");
if(!moondatafile){
cout<<"Cannot open moon data file "<<endl;
exit(-1);
};
for (i=0;i<240000;i++){
moondatafile>>jdtemp>>ra_m[i]>>decl_m[i];
};
moondatafile.close();

ifstream sundatafile("sun_position.dat");
if(!sundatafile){
cout<<"Cannot open sun data file "<<endl;
exit(-1);
};
for (i=0;i<10000;i++){
sundatafile>>jdtemp>>ra_s[i]>>decl_s[i];
};
sundatafile.close();


nevent=data->GetEntries();
cout<<" number of events "<<nevent<<endl;
cout<<fixed;
for (j=0;j<nevent;j++)
//for (j=0;j<1000000;j++)
{data->GetEntry(j);

index=int((jd-2447527.5)*24.);
delx=(jd-(double(index)/24.+2447527.5))/.041667;
temp=ra_m[index+1]-ra_m[index];
if(temp<0.)temp=temp+360.;
ramoon=ra_m[index]+delx*temp;
decmoon=decl_m[index]+delx*(decl_m[index+1]-decl_m[index]);

index=int(jd-2447527.5);
delx=jd-(double(index)+2447527.5);
temp=ra_s[index+1]-ra_s[index];
if(temp<0.)temp=temp+360.;
rasun=ra_s[index]+delx*temp;
decsun=decl_s[index]+delx*(decl_s[index+1]-decl_s[index]);

moonangle=a.calc_angle(ra1,dec1,ramoon,decmoon);

if(moonangle>0)h1->Fill(moonangle,1.0/moonangle);
if(j%1000000==0){
cout<<j<<" "<<jd<<" "<<runevent<<" "<<ra1<<" "<<dec1<<" "<<ramoon<<" "<<decmoon
<<" "<<rasun<<" "<<decsun<<endl;
};
};
h1->Draw("E");

return;

}