#include <fstream>
#include <TMatrixD.h>
#include <TH1D.h>
#include <TMinuit.h>

class clase{

private:

  TObjArray *calibradoIx, *calibradoIy, *linearityIx, *linearityIy;
  //TH1F *residx, *residy, *deflectx, *deflecty;
  char calibrado, linearity, fout;
  double iprxm, iprym, slope, cortex, srpx, srpy, sigmaxtheta, sigmaytheta, residxm, residym, distancia, nfilas, pfilas, fnormresidx, efnormresidx, fnormresidy, efnormresidy, sprx, esprx, spry, espry, sigmaresidx, esigmaresidx, sigmaresidy, esigmaresidy, fnormdeflectx, efnormdeflectx, fnormdeflecty, efnormdeflecty, deflectx, edeflectx, deflecty, edeflecty, sigmadeflectx, esigmadeflectx, sigmadeflecty, esigmadeflecty;
  double calibradox[100], calibradoy[100], fnormcalx[100], ecalibradox[100], ecalibradoy[100], fnormcaly[100], linearityx[100], linearityy[100], elinearityx[100], elinearityy[100], fnormlinx[100], fnormliny[100], sigmacalx[100], sigmacaly[100], sigmalinx[100], sigmaliny[100], plataformax, plataformay, efnormcalx[100], efnormcaly[100], efnormlinx[100], efnormliny[100], esigmacalx[100], esigmacaly[100], esigmalinx[100], esigmaliny[100], residslopex, residslopey, residcortex, residcortey, eresidslopex, eresidslopey, eresidcortex, eresidcortey, residuosx, residuosy, Intensidad[64], sigmax[64], tiras[64];

  enum state {init, full, linok, defok};

public:

  clase();

  void lee(char *name1, char *name2,  Int_t linmax_name1, Int_t linmax_name2);

  void calibradoxpos();

  void calibradoypos();

  void linearityxpos();

  void linearityypos();

  void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);

  Double_t func(Double_t tiras,Double_t par[3]);

};//class

clase::clase(){}

void clase::lee(char *name1, char *name2,  Int_t linmax_name1, Int_t linmax_name2){
  
  char hname1[100];
  char hname2[100];
  char hname3[100];
  char hname4[100];
  
  std::ifstream ino1;
  std::ifstream ino2;
  
  Int_t j=0;
  
  ino1.open(name1);
  ino2.open(name2);
  
  TMatrixD datos(linmax_name1,8);
  TMatrixD centroides(linmax_name2,2);
  TH1D *h1, *h2, *h3, *h4;


  sprintf(hname1,"hola");
  sprintf(hname2,"name2");
  sprintf(hname3,"name3");
  sprintf(hname4,"name4");
  
  
  calibradoIx=new TObjArray(linmax_name2-1,0);
  calibradoIy=new TObjArray(linmax_name2-1,0);
  linearityIx=new TObjArray(linmax_name2-1,0);
  linearityIy=new TObjArray(linmax_name2-1,0);
  
  for(int i=0;i<linmax_name2;i++){
    
    sprintf(hname1,"calIx%i",i);
    sprintf(hname2,"calIy%i",i);
    sprintf(hname3,"linIx%i",i);
    sprintf(hname4,"linIy%i",i);
    h1= new TH1D(hname1,"Perfil de intensidad de carga coordenada x del sensor de calibrado",64,0,63);
    h2= new TH1D(hname2,"Perfil de intensidad de carga coordenada y del sensor de calibrado",64,0,63);
    h3= new TH1D(hname3,"Perfil de intensidad de carga coordenada x del sensor de linearidad",64,0,63);
    h4= new TH1D(hname4,"Perfil de intensidad de carga coordenada y del sensor de linearidad",64,0,63);
    calibradoIx->Add(h1);
    calibradoIy->Add(h2);
    linearityIx->Add(h3);
    linearityIy->Add(h4);
    
  }//for
  
  for (int i=0;i<linmax_name2;i++) {
    
    ino2 >> centroides(i,0) >> centroides(i,1);
    
  }//for
  
  for (int i=0;i<linmax_name1;i++) {
    
    Double_t m;
    ino1 >> m >> datos(i,0) >> datos(i,1) >> datos(i,2) >> datos(i,3) >> datos(i,4) >> datos(i,5) >> datos(i,6) >> datos(i,7);
    
  }//for

  ino1.close();
  ino2.close();
  
  for(int i=0;i<linmax_name2;i++){
    
    Int_t l=1;
    
    for(j=64*(i);j<64*(i+1);j++){
      
      if (datos(j,1)==0.0) datos(j,1)=0.1;
      if (datos(j,3)==0.0) datos(j,1)=0.1;
      if (datos(j,5)==0.0) datos(j,1)=0.1;
      if (datos(j,7)==0.0) datos(j,1)=0.1;
      
      ((TH1D*)calibradoIx->At(i))->AddBinContent(l,datos(j,0)); //numero de entries solo las pone con fill
      ((TH1D*)calibradoIx->At(i))->SetBinError(l,datos(j,1));
      ((TH1D*)calibradoIy->At(i))->AddBinContent(l,datos(j,2));
      ((TH1D*)calibradoIy->At(i))->SetBinError(l,datos(j,3));
      ((TH1D*)linearityIx->At(i))->AddBinContent(l,datos(j,4));
      ((TH1D*)linearityIx->At(i))->SetBinError(l,datos(j,5));
      ((TH1D*)linearityIy->At(i))->AddBinContent(l,datos(j,6));
      ((TH1D*)linearityIy->At(i))->SetBinError(l,datos(j,7));
      l=l+1;
      
    }//for

  }//for

  ((TH1F*)calibradoIx->At(54))->Draw();
  
}//void

void clase::calibradoxpos(){

Double_t para[3]={557.0,27.0,2.5};
Double_t step[3]={0.5,0.5,0.1};

for(int i=0;i<64;i++){

tiras[i]=i;

}//for

for(int j=0;j<linearityIx->GetLast();j++){

for(int i=0;i<64;i++){

Intensidad[i]=((TH1D*)linearityIx->At(j))->GetBinContent(i+1);
sigmax[i]=((TH1D*)linearityIx->At(j))->GetBinError(i+1);

}//for

Int_t ierflg=0;
Double_t arglist[2];
TMinuit *min=new TMinuit(3);
arglist[0]=1;
min->SetFCN(clase::fcn);
min->mnexcm("SET ERR", arglist ,1,ierflg);
min->mnparm(0,"ctex",para[0],step[0],0,0,ierflg); //Faltan los flags
min->mnparm(1,"mediax",para[1],step[1],0,0,ierflg);
min->mnparm(2,"sigmax",para[2],step[2],0,0,ierflg);
arglist[0]=500;
arglist[1]=1.0;
min->mnexcm("MIGRAD",arglist,2,ierflg);
min->mnpout(0,"ctex",fnormlinx[j],efnormlinx[j],0.,0.,0);
min->mnpout(1,"mediax",linearityx[j],elinearityx[j],0.,0.,1);
min->mnpout(2,"sigmax",sigmalinx[j],esigmalinx[j],0.,0.,2);
    
}//for
  
}//void

void clase::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){

Int_t indice=0;
Double_t maximo=0.0;

for(int i=0;i<64;i++){

if(Intensidad[i]<maximo){

maximo=maximo;

}else{

maximo=Intensidad[i];
indice=i;

}//if

}//for

//calculate chisquare
Double_t chisq = 0;
Double_t delta;
for (Int_t i=indice-5;i<indice+5; i++) {

delta  = (Intensidad[i]-(clase::func(tiras[i],par)))/sigmax[i];
chisq += delta*delta;

}//for

f = chisq;

}//void

Double_t clase::func(Double_t tiras,Double_t *par){

Double_t value=par[0]*exp(-0.5*TMath::Power((tiras-par[1])/(par[2]),2));
return value;

}//Float_t