#define Deco_cxx
#include <fstream>
#include <cmath>
#include "Deco.h"
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"

//#define NOPLOT
void Deco::SetBin(int nbin) {_nbin=nbin;} // setta il binning del plot di P
int  Deco::Bin() {return _nbin;}

void Deco::SetZMin(int zmin) {_zmin=zmin;}// decide la Z min per fare il fit
int  Deco::Zmin() {return _zmin;}

void Deco::SetZMax(int zmax) {_zmax=zmax;}// decide la Z MAX per fare il fit
int  Deco::Zmax() {return _zmax;}

void Deco::SetArea(double area) {_area=area;}// (m2) area dell'esperimento
double  Deco::Area() {return _area;}

void Deco::SetTempo(double tempo) {_tempo=tempo;}//(sec) tempo di volo a regime
double  Deco::Tempo() {return _tempo;}

void Deco::SetGeomagE(double gemin) {_gemin=gemin;}// (GeV) cutoff geomagnetico
double  Deco::GeomagE() {return _gemin;}

// (e) risoluzione in carica (metodo per definire membro)
void Deco::SetZRes(double zres) {  Deco::_zres=zres;}
double  Deco::ZRes() {  return _zres; }



void Deco::SetChi2(double chi2) {_chi2=chi2;} // stora il valore del chi2
double  Deco::Chi2() {return _chi2;}

void Deco::SetNDF(int ndf) {_ndf=ndf;}// stora il numero di gradi di liberta'
int  Deco::NDF() {return _ndf;}

void Deco::Loop() // questo e' creato di default da MakeClass, lo tengo per campione
{
//   In a Root session, you can do:
//      Root > .L Deco.C
//      Root > Deco t
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(i);  // read all branches
//by  b_branchname->GetEntry(i); //read only this branch
   if (fChain == 0) return;

   Int_t nentries = Int_t(fChain->GetEntries());

   Int_t nbytes = 0, nb = 0;
   for (Int_t jentry=0; jentry<nentries;jentry++) {
      Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;
   }
}

void Deco::ZDistrib() // questa e' la distribuzione in Z fatta con la curva di claibr. Presto la cancello!
{
   if (fChain == 0) return;

   Int_t nentries = Int_t(fChain->GetEntries());

   char testo[20];



   TH1F *myh1 = new TH1F("myh1","Z nature",_nbin,0,30);
   TH1F *myh2 = new TH1F("myh2","Z integer",_nbin,0,30);
   TH1F *myh3 = new TH1F("myh3","Z+5",_nbin,0,30);
   TH1F *myh4 = new TH1F("myh4","Z integer +5",_nbin,0,30);

   double areamin=100000,areamax=0;
   Int_t nbytes = 0, nb = 0;
   for (Int_t jentry=0; jentry<nentries;jentry++) {
      Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;

      myh1->Fill(tZ);
      myh2->Fill(int(tZ));
      myh3->Fill(tZ+5);
      myh4->Fill(int(tZ)+5);

   }
   TCanvas *nc = new TCanvas("nc");


   
   nc->Divide(2,2);

   nc_1->cd();   myh1->Draw();
   nc_2->cd();   myh2->Draw();
   nc_3->cd();   myh3->Draw();
   nc_4->cd();   myh4->Draw();
}

void Deco::PDistrib(double pmin, double pmax) // grafico di P
{
   if (fChain == 0) return;




   _pmin = pmin; _pmax=pmax;
   Int_t nentries = Int_t(fChain->GetEntries());
   
   TH1F *hp = new TH1F("hp","P distrib",_nbin,pmin,pmax);

   double areamin=100000,areamax=0;
   Int_t nbytes = 0, nb = 0;
   for (Int_t jentry=0; jentry<nentries;jentry++) {
      Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;

      hp->Fill(tp); 
   }

#ifndef NOPLOT
   //   TCanvas *nc = new TCanvas("nc");
   hp->Draw(""); 
#endif
   TFile *hfile = new TFile("Pdist.root","RECREATE","DATA");    
   hp->Write();
   hfile->Close();

}

double Deco::p(int Z)
{
  double Z0 =8.71683e-01;
  double Z1 =3.34011e-02;
  double Z2 =9.80098e-05;
  double Z3 =1.43241e-05;
  
  double zz = Z;
  double pp=  Z0
    +Z1*pow(zz,1)
    +Z2*pow(zz,2)
    +Z3*pow(zz,3);
    
  return pp;
}

double Deco::dp(int Z)
{


  double Z0 =8.71683e-01;
  double Z1 =3.34011e-02;
  double Z2 =9.80098e-05;
  double Z3 =1.43241e-05;
  
  double zz = Z;
  double dzz =_zres;
  double dpp= (
    +Z1
    +2*Z2*pow(zz,1)
    +3*Z3*pow(zz,2)
    )*dzz;
    
  return dpp;
}

TF1 *Deco::gauss(double A, double mean, double sigma)
{
  char ffu[1000];
  sprintf(ffu,"%f*exp(-(pow(x-%f,2))/(2*pow(%f,2)))",A,mean,sigma); 
  TF1 *fu = new TF1("gauss",ffu,1,4.2);
  return fu;
}



double pfunc24(Double_t *x, Double_t *par)
{
  double fitval=0;
  int i=0;
  while(i<24)
    {
  fitval += par[i]*TMath::Exp(-TMath::pow(x[0]-Deco::p(i+5),2)
			      /(2*TMath::pow(Deco::dp(i+5),2)));
  i++;
    }
  return fitval;
}

double pfunc7(Double_t *x, Double_t *par)
{
  double fitval=0;
  int i=2;
  int limite = par[0]+2;
  while(i<limite)
    {
  fitval += par[i]*TMath::Exp(-TMath::pow(x[0]-Deco::p(i+par[1]-2),2)
			      /(2*TMath::pow(Deco::dp(i+par[1]-2),2)));
  i++;
    }
  return fitval;
}

double zfunc7(Double_t *x, Double_t *par)
{
  double fitval=0;
  int i=2;
  int limite = par[0]+2;
  while(i<limite)
    {
  fitval += par[i]*TMath::Exp(-TMath::pow(x[0]-i-par[1]+2,2)
			      /(2*TMath::pow(Deco::ZRes(),2)));
  i++;
    }
  return fitval;
}



void Deco::FitP(int npar)
{

  TFile *f = new TFile("Pdist.root"); // file containing the p plot
  TH1F *data = (TH1F*)f->Get("hp"); 
#ifndef NOPLOT
  data->Draw("");
#endif

  double zmin = double(_zmin)-0.5; // to center Z histo bars
  double zmax = double(_zmax)+0.5;// idem
  int nozbin = _zmax-_zmin+1;// no. of bins
  int limite=nozbin+2;// maximum limit 

  char funame[20];
  sprintf(funame,"pfunc%d",npar);

  TF1 dummy("dummy","gaus",0,1); //keep this
  TF1 *function[25];


  function[7]= new TF1("pfunc7",pfunc7,_pmin,_pmax,limite); // varing fitting function
  function[24]= new TF1("pfunc24",pfunc24,_pmin,_pmax,limite);
  
  function[npar]->SetParameters(limite,60);

  function[npar]->SetParameter(0,nozbin);
  function[npar]->SetParLimits(0,77,77);

  function[npar]->SetParameter(1,_zmin);
  function[npar]->SetParLimits(1,77,77);
  

  data->Fit(funame, "R"); 

  TF1 *fp = data->GetFunction(funame);
//  
//  
//  // plot the gaussians

  int neg_par=0;
  double pa=0;
  double zz=0,acc=0;


  TH1F *hZ = new TH1F("hZ","",nozbin,zmin,zmax);//isto distribuzione di carica
  int i=2;// this start form 2 since the first 2 parameters are fixed.
  while(i<limite)
    {
      pa= fp->GetParameter(i);  
      if(pa<0) neg_par++;
      TF1 *fff = gauss(pa,p(i+_zmin-2),dp(i+_zmin-2));
      fff->SetLineColor(2);
      fff->Draw("same");

      hZ->SetBinContent(i-1,pa);
      i++;
    }
  
  cout<<"Number of Gaussians with Negative amplitude: "<<neg_par<<endl;
  SetChi2(fp->GetChisquare());
  SetNDF(fp->GetNDF());  
//  delete data;


  TFile *hZfile = new TFile("Zdist.root","RECREATE","DATA");    
  data->Write();
  hZ->Write();
  fp->Write();
  hZfile->Close();
}

void Deco::ZData()
{
  TFile *Zf = new TFile("Zdist.root");
  TH1F *hZc = (TH1F*)Zf->Get("hZ");                  
#ifndef NOPLOT
  hZc->SetMarkerStyle(28);
  hZc->SetMarkerSize(0.9);
  hZc->Draw("ap");
#endif
  
}

void Deco::ZTrueDo() // charge-distribution in letteratura
{
  /* crea la distribuzione in Z prendendo i dati dalla letteratura e la pesa
  con l'accettanza del rivelatore. Scrive il risultato nel file ZFdist.root */


  double hmin=_zmin-.5; // estremi del plot consistenti
  double hmax=_zmax+.5;// con il range di Z da fittare 
  int hbin = _zmax-_zmin+1;
  TH1F *hZF = new TH1F("hZF","",hbin,hmin,hmax);//histo distribuzione di carica
  ifstream f_flusso("flussi_relativi_cr.dat");
  ifstream f_accept("accept.dat");
  double zz=0,za=0,fl=0,ac=0;
  int i=1;
  while(f_flusso>>zz>>fl){
    //    cout<<zz<<" "<<_zmin<<endl;
    if(zz>=_zmin && zz<=_zmax){
      while(za != int(zz)){f_accept>>za>>ac;}
      
      int zbin=int(zz)-_zmin+1;
      hZF->SetBinContent(zbin,fl*ac);
    }
  }
  
#ifndef NOPLOT
  //  TCanvas *ZFdis = new TCanvas ("ZFdis");  
  hZF->SetMarkerStyle(29);
  hZF->SetMarkerSize(0.9);
  hZF->Draw("aps");
#endif
  
  TFile *hZFfile = new TFile("ZFdist.root","RECREATE","DATA");    
  hZF->Write();
  hZFfile->Close();

}

void Deco::ZTrue() // recupera e mostra la distribuzione in letteratura
{
  TFile *ZF = new TFile("ZFdist.root");
  TH1F *hZFc = (TH1F*)ZF->Get("hZF");                  
#ifndef NOPLOT
  hZFc->SetMarkerStyle(28);
  hZFc->SetMarkerSize(0.9);
  hZFc->Draw("ap");
#endif
}


void Deco::PTrueConvolve(char* option,int nvb) // Spettro del numero di eventi aspettati
{

  int nozbin = _zmax-_zmin+1;// no. of Zs to consider
  int limite=nozbin+2;// maximum limit of parameter (all the Z + 2 extra info)

  int Zid[100],zsca;
  int is=1;

  double *accettanza[3];
  double acc1[100],acc2[100],acc3[100];
  double fi0[100],gamma[100],nucleoni[100];
  double GEmin=_gemin; // minimum energy (GeV) to avoid the geomagnetic cutoff 

  ifstream spe("spectra.dat");
  ifstream acc("accept3.dat");
  while(spe>>Zid[is]>>nucleoni[is]>>fi0[is]>>gamma[is]){ is++;}
  is=1;
  while(acc>>zsca>>acc1[is]>>acc2[is]>>acc3[is]){ is++;}

  spe.close();
  acc.close();

  accettanza[0]=acc1;
  accettanza[1]=acc2;
  accettanza[2]=acc3;


  TF1 *gausstot= new TF1("pfunc7",pfunc7,_pmin,_pmax,limite); 
  //  gaustot->SetParameters(limite,60);
  gausstot->SetParameter(0,nozbin);
  gausstot->SetParLimits(0,77,77);
  gausstot->SetParameter(1,_zmin);
  gausstot->SetParLimits(1,77,77);

  int iset=0;
  
  while(iset<nozbin){
    
    double esponente =gamma[_zmin+iset]-1; 
    double norm= (
		  fi0[_zmin+iset]
		  *pow(1e3,esponente)
		  *pow(GEmin,-esponente)/(esponente)
		  )
      /sqrt(2*acos(-1)*dp(_zmin+iset))
      /pow(nucleoni[_zmin+iset],esponente)
      *_area
      *_tempo
      *accettanza[nvb][_zmin+iset];
    
    //????      //*TMath::pow(GEmin,gamma[_zmin+iset])
    
    if(_zmin+iset==26)cout<<"ferri "<<norm<<endl;
    
    gausstot->SetParameter(iset+2,norm);
    //    gausstot->SetParLimits(iset+2,77,77);
    
    iset++;
  }
  cout<<"prova1 "<<iset+2<<endl;  


    gausstot->Print();
    gausstot->SetTitle();
    gausstot->Draw(option);
    
   
  // end PTrueConvolve
}

void Deco::ZTrueConvolve(char* option, int nvb) // Spettro del numero di eventi aspettati
{


  int nozbin = _zmax-_zmin+1;// no. of Zs to consider
  int limite=nozbin+2;// maximum limit of parameter (all the Z + 2 extra info)

  int Zid[100],zsca;
  int is=1;

  double *accettanza[3];
  double acc1[100],acc2[100],acc3[100];
  double fi0[100],gamma[100],nucleoni[100];
  double GEmin=_gemin; // minimum energy (GeV) to avoid the geomagnetic cutoff 

  ifstream spe("spectra.dat");
  ifstream acc("accept3.dat");
  while(spe>>Zid[is]>>nucleoni[is]>>fi0[is]>>gamma[is]){ is++;}
  is=1;
  while(acc>>zsca>>acc1[is]>>acc2[is]>>acc3[is]){ is++;}

  spe.close();
  acc.close();

  accettanza[0]=acc1;
  accettanza[1]=acc2;
  accettanza[2]=acc3;

  TF1 *fgaus= new TF1("zfunc7",zfunc7,_zmin,_zmax,limite); 
  //  gaustot->SetParameters(limite,60);
  fgaus->SetParameter(0,nozbin);
  fgaus->SetParLimits(0,77,77);
  fgaus->SetParameter(1,_zmin);
  fgaus->SetParLimits(1,77,77);

  int iset=0;
  
  while(iset<nozbin){
    
    double esponente =gamma[_zmin+iset]-1; 
    double norm= (
		  fi0[_zmin+iset]
		  *pow(1e3,esponente)
		  *pow(GEmin,-esponente)/(esponente)
		  )
      /sqrt(2*acos(-1)*dp(_zmin+iset))
      /pow(nucleoni[_zmin+iset],esponente)
      *_area
      *_tempo
      *accettanza[nvb][_zmin+iset];
    
    //????      //*TMath::pow(GEmin,gamma[_zmin+iset])
    
    if(_zmin+iset==26)cout<<"ferri "<<norm<<endl;
    
    fgaus->SetParameter(iset+2,norm);
    //    fgaus->SetParLimits(iset+2,77,77);
    
    iset++;
  }
  cout<<"prova1 "<<iset+2<<endl;  
  
  
  
    fgaus->Print();
    fgaus->SetTitle();
    fgaus->Draw(option);
    
  //end ZTrueConvolve
}



void Deco::BinCompZ(int elm)
{

  TFile *Zf = new TFile("Zdist.root");
  TFile *ZFf = new TFile("ZFdist.root");


  TH1F *hZc = (TH1F*)Zf->Get("hZ");                  
  TH1F *hZFc = (TH1F*)ZFf->Get("hZF");                  
  
  int NumeroZ=elm; //the element (Z number) whose bin will be used to normalize the histos

  // Normalization of data hZ histo
//  int MaxBin= hZc->GetMaximumBin();
//  double MaximumUp = hZc->GetMaximum();

  int MaxBin = NumeroZ-4;
  double MaximumUp=hZc->GetBinContent(MaxBin);


  TF1* fu = new TF1("fu","1",0,30);

  //  hZc->Divide(fu,MaximumUp);

  // True Flux histo

  double MaxZFc=hZFc->GetBinContent(MaxBin);
  cout<<"Bin + alto: "<<MaxBin<<" Valore: "<<MaxZFc<<endl;
  hZFc->Divide(fu,MaxZFc/MaximumUp);

  
  // -------------------------------------------------------- HISTOS PLOTTING

  TH1F *histo[2];
  histo[0]=hZc;
  histo[1]=hZFc;

  int quale;
  double maximo=0;
  double maxcont=0;
  for(int imax=0;imax<2;imax++)
    {
      maxcont=histo[imax]->GetMaximum();
      if(maxcont>maximo) 
	{
	  maximo=maxcont;
	  quale=imax;
	}
    }
  
//
//
  histo[0]->SetMarkerStyle(20);
  histo[0]->SetMarkerSize(0.9);
  
  histo[1]->SetLineColor(2);
  histo[1]->SetLineWidth(4);
  histo[1]->SetFillColor(2);
  histo[1]->SetFillStyle(3003);


  
  //  ch = new TCanvas("ch","histo");
  
  if(quale==0)
    histo[quale]->Draw("aelp");
  else if (quale==1)
    histo[quale]->Draw("a");
  
  histo[1]->Draw("SAME");
  histo[0]->Draw("elpSAME");
  
  
  
  TLatex *tox[30];
  tox[0]= new TLatex(6,histo[quale]->GetBinContent(7),"C");
  tox[1]= new TLatex(8,histo[quale]->GetBinContent(9),"0");
  tox[2]= new TLatex(14,histo[quale]->GetBinContent(15),"Si");
  tox[3]= new TLatex(26,histo[quale]->GetBinContent(27),"Fe");
  
  for(int ie=0;ie<4;ie++)  
    tox[ie]->Draw();
  

}


void Deco::EvCompZ()
{

  TFile *Zf = new TFile("Zdist.root");
  TFile *ZFf = new TFile("ZFdist.root");


  TH1F *hZc = (TH1F*)Zf->Get("hZ");                  
  TH1F *hZFc = (TH1F*)ZFf->Get("hZF");                  
  
  double NormFact = hZc->Integral()/hZFc->Integral();

  cout<<"Number of Data Events: "<<hZc->Integral()<<endl;

  hZFc->Scale(NormFact);

  
  // -------------------------------------------------------- HISTOS PLOTTING

  TH1F *histo[2];
  histo[0]=hZc;
  histo[1]=hZFc;

  int quale;
  double maximo=0; 
  double minimo=1e6;

  double maxcont=0, mincont=0;
  for(int imax=0;imax<2;imax++)
    {
      maxcont=histo[imax]->GetMaximum();
      mincont=histo[imax]->GetMinimum();
      if(maxcont>maximo) maximo=maxcont;
      if(mincont<minimo) minimo=mincont;
    }

  delete [] histo;

  double hmin=_zmin-.5;
  double hmax=_zmax+.5;
  int hbin = _zmax-_zmin+1;
  TH1F *hcomp = new TH1F("hcomp","",hbin,hmin,hmax);//isto distribuzione di carica
  hcomp->SetMaximum(maximo);
  hcomp->SetMinimum(minimo);

  hZc->SetMarkerStyle(20);
  hZc->SetMarkerSize(0.9);
  
  hZFc->SetLineColor(2);
  hZFc->SetLineWidth(4);
  hZFc->SetFillColor(2);
  hZFc->SetFillStyle(3003);

  
  hcomp->Draw();
  hZFc->Draw("SAME");
  hZc->Draw("elpSAME");

  


  TLatex *tox[30];
  tox[0]= new TLatex(6,hcomp->GetBinContent(7),"C");
  tox[1]= new TLatex(8,hcomp->GetBinContent(9),"0");
  tox[2]= new TLatex(14,hcomp->GetBinContent(15),"Si");
  tox[3]= new TLatex(26,hcomp->GetBinContent(27),"Fe");
  
  for(int ie=0;ie<4;ie++)  
    tox[ie]->Draw();
  

}











//// -------------- Rimasugli obsoleti

//void Deco::EvCompZ()
//{
//
//  TFile *Zf = new TFile("Zdist.root");
//  TFile *ZFf = new TFile("ZFdist.root");
//
//
//  TH1F *hZc = (TH1F*)Zf->Get("hZ");                  
//  TH1F *hZFc = (TH1F*)ZFf->Get("hZF");                  
//  
//  double NormFact = hZc->Integral()/hZFc->Integral();
//
//  cout<<"Number of Data Events: "<<hZc->Integral()<<endl;
//
//  hZFc->Scale(NormFact);
//
//  
//  // -------------------------------------------------------- HISTOS PLOTTING
//
//  TH1F *histo[2];
//  histo[0]=hZc;
//  histo[1]=hZFc;
//
//  int quale;
//  double maximo=0; minimo=1e6;
//
//  double maxcont=0; mincont=0;
//  for(int imax=0;imax<2;imax++)
//    {
//	maxcont=histo[imax]->GetMaximum();
//	mincont=histo[imax]->GetMinimum();
//	if(maxcont>maximo) 
//	  {
//	    maximo=maxcont;
//	    quale=imax;
//	  }
//    }
//  cout<<"Maximum height of histo "<<maxcont<<endl; 
//
//
//
//  histo[0]->SetMarkerStyle(20);
//  histo[0]->SetMarkerSize(0.9);
//  
//  histo[1]->SetLineColor(2);
//  histo[1]->SetLineWidth(4);
//  histo[1]->SetFillColor(2);
//  histo[1]->SetFillStyle(3003);
//
//  
//  //  ch = new TCanvas("ch","histo");
//  
//  if(quale==0)
//    histo[quale]->Draw("aelp");
//  else if (quale==1)
//    histo[quale]->Draw("a");
//
//  histo[1]->Draw("SAME");
//  histo[0]->Draw("elpSAME");
//
//  
//
//
//  TLatex *tox[30];
//  tox[0]= new TLatex(6,histo[quale]->GetBinContent(7),"C");
//  tox[1]= new TLatex(8,histo[quale]->GetBinContent(9),"0");
//  tox[2]= new TLatex(14,histo[quale]->GetBinContent(15),"Si");
//  tox[3]= new TLatex(26,histo[quale]->GetBinContent(27),"Fe");
//  
//  for(int ie=0;ie<4;ie++)  
//    tox[ie]->Draw();
//  
//
//}