#include <string>
#include <fstream>
#include "TRandom.h"
// 
// 

void fit_IndiceDiRifrazionePrisma(string dataFile="SIMUL")
{
    static const Double_t degtorad  = 0.017453293;
    static const Double_t radtodeg  = 57.29577951;

    string fileNameString = "";

    if (dataFile=="")
    {
        std::cout<<"Select an input data file as argument \n Example: \n root[0] .x fit_IndiceDiRifrazionePrisma.C(\"mydata.data\"); \n Otherwise \n root[0] .x fit_IndiceDiRifrazionePrisma.C(\"TEST\") to run this macro on test data"<<std::endl;
        return;
    }
    else if (dataFile == "TEST")
    {
        fileNameString = "/afs/le.infn.it/user/s/spagnolo/html/allow_listing/LabFisicaGenerale/RootMacros/testData/prisma.dat";
    }
    else if (dataFile == "SIMUL")
    {
        simul_Data();
        return;
    }
    else  fileNameString  = dataFile;
    
    
   //
   // To see the output of this macro, click  here. 
   //
   gROOT->Reset();
   //c1 = new TCanvas("c1","Efficiency plateau",200,10,700,500);
   c1 = new TCanvas("c1","Punti Coniugati di una Lente Sottile");


   gPad->SetFillColor(0);
   gPad->SetGrid();
   gPad->GetFrame()->SetFillColor(21);
   gPad->GetFrame()->SetBorderSize(12);

   gStyle->SetOptFit(1112);

   // Reading input data file (x, dx, y, dy)
   char* fileName = fileNameString.c_str();
   cout<<" *** Reading data from file <"<<fileName<<">"<<endl;
   ifstream inputFile;
   inputFile.open(fileName);
   Int_t nP = 0;
   Float_t xVector[100];
   Float_t dxVector[100];
   Float_t yVector[100];
   Float_t dyVector[100];
   while (1) {
       Float_t  xvar = 0; Float_t dxvar = 0.; Float_t yvar = 0; Float_t dyvar = 0.;
       inputFile >> xvar >> yvar >> dxvar >> dyvar;
       if (!inputFile.good()) break;
       cout<<" measurement n. "<<nP+1
           <<"   --- x = "<<xvar
           <<" y = "<<yvar<<" --- errors on x "<<dxvar<<" on y "<<dyvar<<endl;
       // 
       xVector[nP]      = xvar;
       yVector[nP]      = yvar;
       dxVector[nP]     = dxvar;
       dyVector[nP]     = dyvar;
       //
       nP++;
   }
   std::cout<<" *** End of file reached; Number of measurements  = "<<nP<<std::endl;
   if (nP < 1) return;
   
   
   gr = new TGraphErrors(nP,xVector,yVector,dxVector,dyVector);
   string MyTitle = "Indice di Rifrazione del Prisma";
   gr->SetTitle(MyTitle.c_str());
   gr->SetMarkerColor(4);
   gr->SetMarkerStyle(20);
   gr->SetMarkerSize(1);
   gr->GetXaxis()->SetTitle("Angolo di incidenza / deg");
   gr->GetYaxis()->SetTitle("Angolo di deviazione / deg");
   //gr->Draw("ALP");
   gr->Draw("AP");
   cout<<" ************************************************* Plot drawn "<<endl;


   // definiamo una nostra funzione entro il range definito da lowerLimit e upperLimit
   Double_t lowerLimit = 30.0;
   Double_t upperLimit =  90.;
   // la funzione myfunction e` denominate "f1" - ha 1 parametri - e` definita entro i limiti fissati
   TF1 *f1 = new TF1("f1",myfunction,lowerLimit,upperLimit,2);

   // fissiamo i valori (iniziali) dei parametri 
   f1->SetParameter(0,60.);
   f1->SetParameter(1,1.5);
   // fissiamo dei limiti per i valori del parametri 
   //f1->SetParLimits(0,5.,20.);
   f1->SetLineColor(kRed);
   //f1->Draw("L");

   //   return;
   
   
   //fitta entro il range (opzione "R") specificato nella definizione della funzione
   gr->Fit("f1","R"); 
   cout<<" ************************************************* Fit done "<<endl;

   c1->Update();
   
}

void simul_Data()
{
    static const Double_t degtorad  = 0.017453293;
    static const Double_t radtodeg  = 57.29577951;

    std::string fileNameString = "simulation_prisma.dat";
    char* fileName = fileNameString.c_str();
    cout<<" *** simulating data in file <"<<fileName<<">"<<endl;
    std::ofstream outputFile(fileName);
    //outputFile.open(fileName);

    Int_t nP = 50;
    Double_t theta_start = 30.;
    Double_t theta_stop = 80.;
    Double_t step = (theta_stop-theta_start)/Double_t(nP-1);
    
    Double_t theta = theta_start;
    Double_t x[1];
    Double_t par[2];
    par[0] = 60.;
    par[1] = 1.5;
    Double_t sigma_x = 0.03;
    Double_t sigma_y = 0.05;
    Double_t measured_sigma_x = sigma_x;
    Double_t measured_sigma_y = sigma_y;
    Double_t nsig_x = 0.;
    Double_t nsig_y = 0.;
    for (int i=0; i<nP; i++)
    {
        x[0] = theta;
        Double_t delta = myfunction(x, par);
        //        std::cout<<i<<" --  theta - delta  "<<x[0]<< "          "<<delta<<std::endl;
        
        nsig_x = gRandom->Gaus();
        nsig_y = gRandom->Gaus();

        delta = myfunction(x, par);
        
        x[0] = x[0]   + nsig_x * sigma_x;
        delta = delta + nsig_y * sigma_y;
        std::cout<<i<<" --  after smearing "<<x[0]<<"          "<<delta<<" n_sigma_x,y "<<nsig_x<<" "<<nsig_y<<std::endl;
        outputFile << " "<<x[0]<<" " << delta <<" "<< measured_sigma_x <<" "<< measured_sigma_y <<std::endl;
        
        theta = theta+step;
    }
    outputFile.close();

    return;
    
}

Double_t myfunction(Double_t *x, Double_t *par)
{
    const Double_t degtorad  = 0.017453293;
    const Double_t radtodeg  = 57.29577951;
    //     par[0] = angolo di rifrangenza (deg)
    //     par[1] = indice di rifrazione

    //     x[0] = angolo di rifrangenza (deg)
    Double_t theta = x[0]*degtorad; // in radiants here 
    
    Double_t A = par[0];
    A = degtorad*A; // convert to radians

    Double_t n = par[1];
    //std::cout<<" x, A, theta, n "<<x[0]<<" in rad = "<<theta<<" "<<A<<" "<<n<<std::endl;

    
    Double_t delta = theta - A + asin( n*sin( A - asin( sin(theta)/n ) ) );
    delta = radtodeg * delta; // in deg
    //std::cout<<"delta = "<<delta<<std::endl;

    //std::cout<<" "<<x[0]<<"      "<<delta<<std::endl;

    
    return delta;
}