// system includes

#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <iostream>
#include <fstream>

using namespace std;

// root includes

#include "TApplication.h"
#include "TLine.h"
#include <TROOT.h>
#include "TStyle.h"
#include "TF2.h"
#include "TF1.h"
#include "TH2F.h"
#include "TPaveLabel.h"
#include "TCanvas.h"
#include "TPostScript.h"
#include "TMultiGraph.h"
#include "TGraph.h"
#include "TText.h"
#include "TArrow.h"
#include "THStack.h"

//#include "TH3F.h" //EXPERIMENTAR 3D

#include "TMath.h" //included because compile error 20100313

// user includes

#include "userin.h"
#include "TCurves.h"


// user defined expressions

# define sq(x) (x)*(x) // x squared

#define GAMMA    22
#define ELECTRON 11
#define POSITRON -11
#define PROTON   2212
#define NEUTRON  2112
#define C12      1000060120
#define He3      1000020030
#define He4      1000020040
#define T        1000010030
#define D        1000010020
#define PIplus   211
#define PIminus  -211
#define PI0      111

#define TWOP_RESTMASS 1876.6 // MeV
#define nRMsq 882848.16      // MeV^2
#define nRM   939.6          // MeV

#define PI 3.14159265
#define CONST 4*PI*300*300

//Added 100502 MC @HOME
#define PROTON_RESTMASS 1.67262158e-27 //Kg

 
// user defined global variables

TCanvas *gC;     // canvas 
//TCanvas *gT;     // timing canvas 
TGraph  *gr1;    // graph
TGraph  *gr2;    // graph
TGraph  *grav1;  // graph average
TGraph  *grav2;  // graph average
TGraph  *gren1;  // graph energy
TGraph  *gren2;  // graph energy
TGraph  *grtm;   // graph time
TH1F    *h1tm;   // histogram with coincidence time
TH1F    *h1tof;   // histogram with coincidence time
TH1F    *h1NeutronEnergy; // histogram with energy (qdc)
TH1F    *h1LatScat; // histogram with energy (qdc)
TH2F    *h2FinalGamma; // histogram with energy (qdc)
TH1F    *h1FinalGamma; // histogram with energy (qdc)
TH2F    *h2Stop2D; // histogram with energy (qdc)
TH1F    *h1en_1t;// histogram with pulse height (trigger window)
TH1F    *h1en_1tbck;// histogram with pulse height (bck window)
TH1F    *h1en_2t;// histogram with pulse height (trigger window)
TH1F    *h1en_2tbck;// histogram with pulse height (bck window)
TH1F    *h1en_1slope;// histogram with slope (trigger window)
TH1F    *h1en_2slope;// histogram with slope (trigger window)
TH1F    *h1en_1slope_notrigg;// histogram with slope (trigger window)
TH1F    *h1en_2slope_notrigg;// histogram with slope (trigger window)
TH1F    *h1en_1negslope;// histogram with slope (trigger window)
TH1F    *h1en_1negslope_notrigg;// histogram with slope (trigger window)
THStack *hs;//("hs","test stacked histograms");

TF1 *vfit;

// function definition

Double_t fitgauss(Double_t *x, Double_t *par)
{
   Double_t arg = 0;
   if (par[2] != 0) arg = (x[0] - par[1])/par[2];

   // more parameters, PC070125
   // ori Double_t fitval = par[0]*TMath::Exp(-0.5*arg*arg)+100;
   Double_t fitval = par[0]*TMath::Exp(-0.5*arg*arg)+par[3]+par[4]*x[0];

   return fitval;
}


Double_t fitLatScat(Double_t *x, Double_t *par)
{
 
  Double_t VamosVoigt = TMath::Voigt(x[0],par[0], par[1], 4);
  
  return VamosVoigt;

}

int main(int argc, char **argv)
{
  Int_t       i_t, j_t,          // loop variables
               nNeutrons,         //number of produced neutrons
              //tot_i_t=109,    // total iterations on final events
              tot_i_t=972252,    // GAMMAs total iterations on final events
              //tot_i_t=12500,    // total iterations on final events
              //tot_i_t=5000,    // total iterations on final events
              *partType;        // particle type (HEP encoding)
  char        lmdfile[100],     // waveform file to be read
              aux_text[100];    // auxiliary text
  Double_t    *E_i,             // initial energy (MeV)
    //        gammaE,           // energy of final gamma ray (MeV)
    //        gammaP,           // momentum of final gamma ray (MeV)
              KE,               // kinetic energy
              R, Theta,    //Spherical coordinates
              neutronDirection[3],
              Zangle, i_t_RAD,          // polar angle
              //x_,y_,z_,
              *partProp,        // particle properties (r, p)
              *tof;
  FILE        *fh;              // file handle

  int      *profile,
           *profile_proj,
           *dist2D,
           zint, //mu
           xint_10mu, zint_10mu, //10 mu
           *initE,
           *escapeEall,
           *escapeEaccepted,
           *flightTimePart1,
           *flightTimePart2,
           *rest2D;

  double   Phi, area, area_degree, Phi_before_rad, Phi_middle, Phi_middle_rad, Phi_after_rad;
  float    Phi_before;
  int       result;


 printf("Sizeof double=%d, char=%d\n",sizeof(double),sizeof(char));

 // variable initialization and memory allocation

 // profile per mm
 if ( (profile = (int *) calloc(301,sizeof(int))) == NULL ){
   puts("Problem allocating memory for \"profile\". Exiting.\n");
   exit(1);
   }
 // profile per mu
 if ( (profile_proj = (int *) calloc(5001,sizeof(int))) == NULL ){
   puts("Problem allocating memory for \"profile_proj\". Exiting.\n");
   exit(1);
   }
 // 2D distribution per mu
 if ( (dist2D = (int *) calloc(501*301,sizeof(int))) == NULL ){
   puts("Problem allocating memory for \"dist2D\". Exiting.\n");
   exit(1);
   }
 // initE per 100 keV
 if ( (initE = (int *) calloc(201,sizeof(int))) == NULL ){
   puts("Problem allocating memory for \"profile\". Exiting.\n");
   exit(1);
   }
 // escapeE per 10 keV (total)
 if ( (escapeEall = (int *) calloc(100000,sizeof(int))) == NULL ){
   puts("Problem allocating memory for \"profile\". Exiting.\n");
   exit(1);
   }
 // escapeE per 10 keV (accepted)
 if ( (escapeEaccepted = (int *) calloc(100000,sizeof(int))) == NULL ){
   puts("Problem allocating memory for \"profile\". Exiting.\n");
   exit(1);
   }
 // flight time per 100 ps
 if ( (flightTimePart1 = (int *) calloc(500,sizeof(int))) == NULL ){
   puts("Problem allocating memory for \"profile\". Exiting.\n");
   exit(1);
   }
 // flight time per 100 ps
 if ( (flightTimePart2 = (int *) calloc(500,sizeof(int))) == NULL ){
   puts("Problem allocating memory for \"profile\". Exiting.\n");
   exit(1);
   }
 // rest2D per mm
 if ( (rest2D = (int *) calloc(90601,sizeof(int))) == NULL ){
   puts("Problem allocating memory for \"profile\". Exiting.\n");
   exit(1);
   }



 // particle type
 if ( (partType = (int *) calloc(1,sizeof(int))) == NULL ){
   puts("Problem allocating memory for \"partType\". Exiting.\n");
   exit(1);
   }
 // initial energy
 if ( (E_i = (double *) calloc(1,sizeof(double))) == NULL ){
   puts("Problem allocating memory for \"partType\". Exiting.\n");
   exit(1);
   }
 // particle properties
 if ( (partProp = (double *) calloc(6,sizeof(double))) == NULL ){
   puts("Problem allocating memory for \"partProp\". Exiting.\n");
   exit(1);
   }
 // tof
 if ( (tof = (double *) calloc(1,sizeof(double))) == NULL ){
   puts("Problem allocating memory for \"tof\". Exiting.\n");
   exit(1);
   }


 // root handling

 TApplication theApp("App", &argc, argv);
 gROOT->Reset();
 gStyle->SetOptStat(1111111);  //adjust statistics box
 //gStyle->SetOptStat(10);  //adjust statistics box
 gStyle->SetPalette(1);

 // event analysis
 gC = new TCanvas("gC","Neutron production",300,50,1100,900);
 gStyle->SetCanvasColor(0);
 gStyle->SetFrameFillColor(0);
 gC->SetFillColor(0);
 gC->SetGrid();
 gC->Divide(2,2,.0001,.0001);

 ////////////////////////
 // root physics handling
 ////////////////////////

 // histogram with initial particle energy
 h1NeutronEnergy = new TH1F("h1NeutronEnergy","Initial E",2000, 0,20);

 h1LatScat = new TH1F("h1LatScat","Lateral Scattering",1800,0.,180.);

 h1FinalGamma = new TH1F("h1FinalGE","Final Gamma",1800,0.,1.);

 h2Stop2D = new TH2F("h1FinalGE","Final Gamma",7001,-350.,350.,7001,-350.,350.);


 ////////////
 // lmd file
 ////////////

 sprintf(lmdfile,"/home/micaela/g4work/pv100/output/protonWater_18.00MeV_2M_102.lmd");

 cout << "Opening lmd file \"" << lmdfile << "\"."<< endl << flush;

 // open lmd file
 if ( (fh = fopen (lmdfile,"rb")) == NULL ) {
   printf ("\nCannot open file \"%s\".\n", lmdfile);
   printf ("Exiting.\n"); 
   exit(1);
   }

 rewind(fh);

 i_t = 0;

 while (!feof(fh)) {

   // read particle type
   if (fread( partType, sizeof(int), 1, fh) != 1) {
     printf("\n Error reading particle type from file \"%s\"!!!\n",lmdfile);
   }
   // read initial energy
   if (fread( E_i, sizeof(double), 1, fh) != 1) {
     printf("\n Error reading inital energy from file \"%s\"!!!\n",lmdfile);
   }
   // read particle properties
   if (fread( partProp, sizeof(double), 6, fh) != 6) {
     printf("\n Error reading particle properties from file \"%s\"!!!\n",lmdfile);
   }

   if (*partType != -1 && *E_i != 0.) {

     KE = sqrt(sq(*(partProp+3)) + sq(*(partProp+4)) + sq(*(partProp+5)) +
	   nRMsq) - nRM;

     if (*partType==PROTON && KE != 0.){

       zint = (int)floor((*(partProp+2) * 1.e3))+2500; //mu
       xint_10mu = (int)floor((*(partProp+0) * 1.e2))+150; //10mu
       zint_10mu = (int)floor((*(partProp+0) * 1.e2))+250; //10mu
       if (zint >= 0 && zint <= 5000) {
         *(profile_proj+zint) = *(profile_proj+zint)+1;
         if (zint_10mu >= 0 && zint_10mu <= 500 &&
             xint_10mu >= 0 && xint_10mu <= 300 ) {
           *(dist2D+zint_10mu*300+xint_10mu) = *(dist2D+zint_10mu*300+xint_10mu)+1;
	 }         
       }

       R=300.;


       if (*(partProp+2) <= 300){

	 Phi = acos((*(partProp+2))/R)*(180./PI);
	 Phi_before = int(Phi*10.)/10.;              
	 Phi_before_rad = Phi_before*(PI/180.);

	 Phi_middle = (Phi_before+0.05);
	 Phi_after_rad = (Phi_before+0.1)*(PI/180.);   

	 area = 2*PI*R*R*((cos(Phi_before_rad))-(cos(Phi_after_rad)));

	 h1NeutronEnergy->Fill(KE,1.);
	 if (*(partProp+0)>=-2. && *(partProp+0)<2)
	   h2Stop2D->Fill(*(partProp+2),*(partProp+1));
	 h1FinalGamma->Fill(((Phi_middle/180)),(CONST/area));
	 h1LatScat->Fill(Phi_middle,(1/area));
       }	 
     }

   i_t++;

   } // end loop events and read file

 // close input file
 fclose (fh);


 ////////
 // plot
 ////////
 
 // initial energy
 gC->cd(1);
 h1NeutronEnergy->SetXTitle("Energy (MeV)");
 h1NeutronEnergy->SetYTitle("Counts per 100 KeV");
 h1NeutronEnergy->GetXaxis()->CenterTitle();  //align axis title in the center  
 h1NeutronEnergy->GetYaxis()->CenterTitle();
 // //h1initE->SetTitleSize(.7,"Y");  //adjust the size of the axis title
 h1NeutronEnergy->Draw(); // scatter plot


 // detected energy
 gC->cd(2);
 TF1 *fitFcn = new TF1("fitFcn",fitgauss,0.,180.,3);

 h1LatScat->SetXTitle("degrees");
 h1LatScat->SetYTitle("Counts (dN/dTheta)");
 h1LatScat->GetXaxis()->CenterTitle();  //align axis title in the center  
 h1LatScat->GetYaxis()->CenterTitle();
 //h1initE->SetTitleSize(.7,"Y");  //adjust the size of the axis title
 h1LatScat->Draw(); // scatter plot

 
 vfit = new TF1("vfit", "Voigt", 0., 180.);

 h1LatScat->Fit(vfit, "RS");
 vfit->SetLineColor(3);
 vfit->Draw("same");
 

 // detected position and energy
 gC->cd(3);
 h1FinalGamma->SetXTitle("degrees");
 h1FinalGamma->SetYTitle("Counts dN/d(sr)");
 h1FinalGamma->GetXaxis()->CenterTitle();  //align axis title in the center  
 h1FinalGamma->GetYaxis()->CenterTitle();
 h1FinalGamma->Draw();
 
 // detected stopping position 
 gC->cd(4);
 h2Stop2D->SetXTitle("Z (mm)");
 h2Stop2D->SetYTitle("Y (mm)");
 h2Stop2D->SetZTitle("Counts");
 h2Stop2D->GetXaxis()->CenterTitle();  //align axis title in the center  
 h2Stop2D->GetYaxis()->CenterTitle();
 h2Stop2D->GetZaxis()->CenterTitle();
 h2Stop2D->Draw("colZ"); // lego plot

 theApp.Run();

 return 0;
}