#include <string>

// 
// 

void fit_LenteSottile(string dataFile="")
{
    string fileNameString = "";

    if (dataFile=="")
    {
        std::cout<<"Select an input data file as argument \n Example: \n root[0] .x fit_LenteSottile.C(\"mydata.data\"); \n Otherwise \n root[0] .x fit_LenteSottile.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/lente.dat";
    }
    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 = "Lente Sottile: ramo XXXX";
   gr->SetTitle(MyTitle.c_str());
   gr->SetMarkerColor(4);
   gr->SetMarkerStyle(21);
   gr->GetXaxis()->SetTitle("Distanza oggetto / units");
   gr->GetYaxis()->SetTitle("Distanza immagine / units");
   //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 = -20.0;
   Double_t upperLimit =  50.;
   // la funzione myfunction e` denominate "f1" - ha 1 parametri - e` definita entro i limiti fissati
   TF1 *f1 = new TF1("f1",myfunction,lowerLimit,upperLimit,1);

   // fissiamo i valori (iniziali) dei parametri 
   f1->SetParameter(0,9.);
   // fissiamo dei limiti per i valori del parametri 
   //f1->SetParLimits(0,5.,20.);
   f1->SetLineColor(kRed);
   
   //fitta entro il range (opzione "R") specificato nella definizione della funzione
   gr->Fit("f1","R"); 
   cout<<" ************************************************* Fit done "<<endl;

   c1->Update();
}

Double_t myfunction(Double_t *x, Double_t *par)
{
    //     par(0) -> distanza focale
    //     x[0] = distanza oggetto 

    //Double_t distanzaImmagine = 1./(1./par[0] - 1./x[0]);
    Double_t distanzaImmagine = x[0]*par[0]/(x[0] - par[0]);
    
    return distanzaImmagine;
}