// C/C++ headers
//
#include <fstream>
#include <string.h>
#include <iostream>
#include <stdio.h>
#include <ctype.h>
//
// ROOT headers
//

#include <TTree.h>

#include <TSystem.h>
#include <TStopwatch.h>
#include <TSystemDirectory.h>
#include <TGraphAsymmErrors.h>
#include <TString.h>
#include <TFile.h>
#include <TClass.h>
#include <TGraph.h>
#include <TLine.h>
#include <TH1F.h>
#include <TMath.h>
#include <TRandom3.h>
#include <TF1.h>
//
// RunInfo header
//
#include <RunInfo.h>

#include <PamLevel2.h>

#include <TClonesArray.h>
#include "PrimaryInfo.h"

using namespace std;
using TMath::Pi;
//This program generate trajectories, which are going in acceptance


void Usage ();



Bool_t IsInside(Trajectory *t, Double_t fiducial,
		float *xmin, float *xmax,  float *ymin,  float *ymax){

  for(Int_t i = 0; i<t->npoint; i++){
    if( !(((xmin[i]+fiducial)<t->x[i])*(t->x[i]<(xmax[i]-fiducial))) ) return kFALSE;
    if( !(((ymin[i]+fiducial)<t->y[i])*(t->y[i]<(ymax[i]-fiducial))) ) return kFALSE;
    //cout<<"XMIN:"<<xmin[i]
    //<<" XMAX:"<<xmax[i]
    //<<" YMIN:"<<ymin[i]
    //<<" YMAX:"<<ymax[i]
	 //<<" Z:"<<TrkParams::zGF[i]<<endl;
  }
  //t->Dump();
  
  return kTRUE;
}



void DoHistogramFromGraph(TGraph* g, TH1F* htemplate, Double_t factor){
  //TH1F* hout = (TH1F*) htemplate->Clone(name);
  htemplate->Reset();
  int n = htemplate->GetNbinsX();
  for(int bin=1; bin<=n; bin++) {
    htemplate->SetBinContent(bin, g->Eval(htemplate->GetBinCenter(bin))*factor);
    Double_t *gx = g->GetX();
    int p = 0;
    Double_t x = htemplate->GetBinCenter(bin);
    for(p=1; p<g->GetN(); p++)
      if((gx[p-1] <= x) && (gx[p] > x)) break;
    Double_t gey = g->GetErrorY(p-1)*factor;
    //cout<<"gx: "<< hout->GetBinCenter(bin)<<"gy: "<<hout->GetBinContent(bin)<<endl;
    htemplate->SetBinError(bin, gey);
    //cout<<"gx: "<< htemplate->GetBinCenter(bin)<<"gy: "<<htemplate->GetBinContent(bin)<<endl;
  }
}

void BinLogX(TH1*h){
  //
  // Modify histogram h to create an histogram logarithmically binned
  // along the X axis.  It works on normally defined, not
  // filled,histograms.
  //
  TAxis *axis = h->GetXaxis();
  int bins = axis->GetNbins();

  Axis_t xmin = axis->GetXmin();
  Axis_t xmax = axis->GetXmax();

  if (xmin==0)
    {
      xmin = 1e-6;
      cout<<"Warning: low edge of histogram set from 0 to 1e-6"<<endl;
    }
  Axis_t from = TMath::Log10(xmin);
  Axis_t to = TMath::Log10(xmax);

  //   Axis_t width = (xmax - xmin) / bins;
  Axis_t width = (to - from) / bins;
  Axis_t *new_bins = new Axis_t[bins + 1];

  for (int i = 0; i <= bins; i++) {
    new_bins[i] = TMath::Power(10, from + i * width);
  }
  axis->Set(bins, new_bins);
  delete new_bins;
}


int main(int argc, char * argv[] ){
 TStopwatch timer;
 timer.Start();
 
 TString DataPATH  = gSystem->Getenv("PWD");


  TString inputfilename, outputfilename;
  TFile* rfin = 0;
  TFile* rfout = 0;
  TTree* primaries = 0;
  PrimaryInfo *p = 0; 
  Double_t fiducial, Z, Rmin, Rmax, R;
  Int_t PDG;
  fiducial=Z=Rmin=Rmax=R=0.;
  Int_t nev;
  TF1 *frigflat = 0;
  TH1F * flux = new TH1F("flux", "flux", 2000, 0.4, 1000.);

  enum mode {SPECTRA, FLATBIN, SINGLEPOINT};
  mode md;

  TString innf(argv[5]);
  switch(argc){
  case 7:
    if(innf.Contains(".root")) md = SPECTRA; else md = SINGLEPOINT;
    break;
  case 8:
    md = FLATBIN;
    break;
  default:
    Usage(); 
    return 0;
    break;
  }

  nev = atoi(argv[1]);
  Z = atof(argv[2]);
  PDG = atoi(argv[3]);
  fiducial = atof(argv[4]);


  TSystemDirectory *workdir = 0;

  switch(md){
  case SPECTRA:
    inputfilename = TString(argv[5]);
    outputfilename = TString(argv[6]);
    workdir=new TSystemDirectory("workdir",DataPATH);
    rfin = new TFile(inputfilename);
    if (rfin){
      cout<<"Generating spectra from file: "<<inputfilename<<endl;
      //Taking RIGIDITY p.d.f.
     
      // TList * outlist = (TList*)rfin->Get("outlist");
      //flux = (TH1F*)outlist->At(5); //CHANGE THIS IF NEED
      TGraphAsymmErrors *gr1 = (TGraphAsymmErrors*)rfin->Get("gflux");
      //TH1F * temp = new TH1F("template", "template", 40, 0.4, 250.);
      BinLogX(flux);
      DoHistogramFromGraph(gr1, flux, 100000.);
      rfin->Close();
      //cout<<"test"<<flux->GetBinCenter(10)<<endl;
    }
    break; 
  case FLATBIN:
    Rmin = atof(argv[5]);
    Rmax = atof(argv[6]);
    if (Rmin == Rmax) cout<<"ERROR: Rmin = Rmax"<<endl;
    frigflat = new TF1("fmom","1",Rmin,Rmax);
    outputfilename = TString(argv[7]);
    cout<<"Generating flat spectra, Rmin="<<Rmin<<" GV,"<<" Rmax="<<Rmax<<" GV"<<endl;
    break;
  case SINGLEPOINT:
    R = atof(argv[5]);
    outputfilename = TString(argv[6]);
    cout<<"Generating fixed rigidity, R="<<R<<" GV"<<endl;
    break;
  default:
    break;
  } 
 
  //Creating output root-file
  rfout = new TFile(outputfilename.Data(),"recreate");
  primaries = new TTree("primaries","primaries");
  TTree::SetMaxTreeSize(1000*Long64_t(2000000000));
  p = new PrimaryInfo();
  if(!primaries->GetBranch("PRIM"))
    primaries ->Branch("PRIM","PrimaryInfo", &p,32000,99);
  else
    primaries->GetBranch("PRIM")->SetAddress(&p);
  
 


  //Loading field, Initialization
  
  TRandom3 *ran = new TRandom3(0);
  TF1 *ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.);
  ftheta->SetNpx(99999);
  
  cout<<"Load Magnetic Field... "<<endl;
  TrkTrack *track = new TrkTrack();
  track->LoadField(getenv("PAM_FIELD"));
  track->SetTrackingMode(0);
  cout<<"Magnetic Field Loaded. "<<endl;

  //Start

  cout<<"Nevents: "<<nev<<", FIDUCIAL: "<<fiducial<<" (cm), Z:"<<Z<<", PDG:"<<PDG<<endl;   
   
#define ZCAL   13.05  //cm
#define CALS3  10.639 //cm
#define ZSPEBP 2.97   //cm
#define HORMG  22.57  //cm


#define XMINGEN -30.
#define XMAXGEN  30.
#define YMINGEN -30.
#define YMAXGEN  30.
#define ZINP    110.

   setprecision(8);

   float *xmin = TrkParams::xGF_min;
   float *xmax = TrkParams::xGF_max;
   float *ymin = TrkParams::yGF_min;
   float *ymax = TrkParams::yGF_max;
   
   Float_t zin[TrkParams::nGF];
   for (Int_t igf=0; igf<TrkParams::nGF; igf++) zin[igf]=TrkParams::zGF[igf];
   Trajectory* t = new Trajectory(TrkParams::nGF,zin);
   Float_t al[5];

   //Main LOOP
   
   Int_t goodev = 0;
   Int_t totev = 0;
   Double_t mom = 0.;
   Double_t x,y,z,theta,phi;
   while (goodev<nev){
     x = ran->Uniform(XMINGEN, XMAXGEN);
     y = ran->Uniform(YMINGEN, YMAXGEN);
     z = ZINP-ZCAL-CALS3-ZSPEBP-HORMG;
     theta = ftheta->GetRandom();
     phi = ran->Uniform(0.,2.*acos(-1.));
     
     switch(md){
     case SPECTRA:
       mom = Z*flux->GetRandom();
       break;
     case FLATBIN:
       mom = Z*frigflat->GetRandom();
       break;
     case SINGLEPOINT:
       mom = Z*R;
       break;
     default:
       break;
     }
     
     al[0] = x;
     al[1] = y;
     al[2] = sin(theta);
     al[3] = phi;
     al[4] = Z/mom;

     t->DoTrack(al, z);

     if(IsInside(t,fiducial, xmin, xmax, ymin, ymax)){
       p->X0 = x;
       p->Y0 = y;
       p->Z0 = ZINP;
       p->Theta = theta;
       p->Phi = phi;
       p->P0 = mom;
       p->PDG=PDG; //He-4
       primaries->Fill();
       goodev++;
     }
     totev++;
     if( ((totev%10000)==0) || (goodev==nev) ) cout<<"-> "<<totev<<" "<<goodev<<endl;
   }
   primaries->Write();
   if(md==SPECTRA) flux->Write(); 
   rfout->Close();
   
   Double_t gf = Double_t(goodev)/Double_t(totev)*(XMAXGEN-XMINGEN)*(YMAXGEN-YMINGEN)*TMath::Pi()/2.;
   Double_t err_gf = (Double_t)1./TMath::Sqrt(goodev)*gf;
   cout<<"GF="<<gf<<" +- "<<err_gf<<endl;
   
   //delete primaries;
   delete p;
   delete ran;
   delete ftheta;
   delete track;
   delete t;
   return 0;   
}

void Usage(){
  cout<<"NAME"<<endl;
  cout<<"       sp_sen - spectra generator"<<endl;
  cout<<"SYNOPSIS"<<endl;
  cout<<"       sp_gen  NUMPAR Z PDG FIDUCIAL(cm) RIGIDITY(GV) OUTFILE"<<endl; 
  cout<<"       sp_gen  NUMPAR Z PDG FIDUCIAL(cm) RIGMIN(GV) RIGMAX(GV) OUTFILE"<<endl; 
  cout<<"       sp_gen  NUMPAR Z PDG FIDUCIAL(cm) INPUTFILE OUTFILE"<<endl; 
  cout<<"DISCRIPTION"<<endl;
  cout<<"       Generate tracks above the apparatus "<<endl; 
  cout<<"       isotropically and save initial kinematic parameters (X0,Y0,Z0=110., P, Theta, Phi)"<<endl;
  cout<<"       good tracks in output root-file. After this, file could be used to do simulations."<<endl;
  cout<<"       The condition to accept tracks is satisfaction to 14 acceptance planes, defined by"<<endl;
  cout<<"       TrkParams in Level2. This tool is based on DoTrack method of Level2. "<<endl;
  cout<<"       To save data in root-file PrimaryInfo class was developed and compiled in library. "<<endl;
  cout<<"       If later you wish to develop your own script to work with data stored threre, "<<endl;
  cout<<"       you need just load libPrimaryInfo.so shared library"<<endl; 
  cout<<"NOTE"<<endl;
  cout<<"       Use http://pdg.lbl.gov/2007/reviews/montecarlorpp.pdf document to define your"<<endl;
  cout<<"       particle according PDG numeration scheme."<<endl;
}