//////////////////////////////////////////////////////////////////////////////////////////////////
///
/// J. Diaz, Feb-2018
///
/// Program to calibrate a DSSSD in a raw: 16 front and 16 back strips
///
/// Specify input file, detector name and energy values of the calibration sources
///
/// Output as .txt file with the slope and offset of the 16+16 strips calibration
///
/// To execute, in the same directory with the input file:
///	root -l
///	.L MyAutoFit.C
///	MyAutoFit()
///
//////////////////////////////////////////////////////////////////////////////////////////////////

#include "TMath.h"
#include <TChain.h>
#include <TH1F.h>
#include <fstream>
#include <iostream>
#include <cstdlib>
class TF1;
class TH1;
class TAxis;
class TMethodCall;
using namespace std;

int MyAutoFit()
{
    gROOT->Reset();

    //declaration of variables
    int binmax=0;
    double x=0., xmax=0., xmin=0., y=0., aux=0., aux2=0., Mean[4], Sigma[4], FWHM[4];
    double slope, offset, sumx, sumy, sumxx, sumxy, dynrange;
    double res[4];
    char strip[35];
    TH1F *h1 = new TH1F("h1","h1",400.,100.,500.); // specify binning and range of the histogram to plot (where the 4 peaks appear)
    TH1F *h2 = new TH1F("h2","h2",400.,100.,500.); // auxiliar histogram

    char* E[2] = {"A_F_E>>h1","A_B_E>>h1"}; // specify here detector to calibrate: A_F_E, A_B_E, B_F_E, B_B_E...
    double energy[4]={3182.69,5156.59,5486,5804.96}; // energies (keV) of the 4 peaks in the alpha source

    TChain ch = TChain("h101");
    ch.Add("6_00*.root"); // specify the files you're using to calibrate, run 6 is a good one to calibrate DeltaE detectors
				//using the whole run might take a time to execute...
    
    //output file name
    ofstream file("Calout_A.txt");
    
	//loop for front-back strips
	for(int k=0; k<=1; k++){
		
	    if(k==0){file << "#Front"<<'\n';}
	    if(k==1){file << "#Back"<<'\n';}

	    //loop for 16 strips	    
	    for(int i=1; i<=16;i++)
	    {
	        if(k==0){sprintf(strip,"A_FI==%d",i);} // specify here detector name 
	        if(k==1){sprintf(strip,"A_BI==%d",i);} // specify here detector name
		h101->Draw(E[k],strip);

		cout << "*****Strip # " << i << endl;
		
		//loop for 4 peaks
		for(int m=0; m<4; m++){
		    Mean[m] = 0.; Sigma[m] = 0.; FWHM[m] = 0.;
		    
		    binmax = h1->GetMaximumBin();
		    x = h1->GetXaxis()->GetBinCenter(binmax);
		    xmin = x-10; xmax = x+10; // specify here region to fit
		    TF1 *g = new TF1("g","gaus",xmin,xmax);
		    h1->Fit(g,"RQ"); 
		    //Q: quiet mode doesn't show parameters in terminal
		    Mean[m] = g->GetParameter(1);
		    Sigma[m] = g->GetParameter(2);
		    FWHM[m] = 2.355 * Sigma[m];
		    
		    //cout << "========================================" << endl;
		
		    for (int j=0; j<2000; j++){
		        h2->SetBinContent(j,0.);}

		    for (int j=(binmax-10); j<(binmax+10); j++){
		        y=h1->GetBinContent(j);
		        h2->SetBinContent(j,y);}

		    h1->Add(h2,-1);
		} 
		
		//sorting peaks from minor to major
		for(int l=0; l<4; l++){
			for(int n=0; n<4; n++){
				if(Mean[l] < Mean[n]){
					aux=Mean[l]; aux2=FWHM[l];
					Mean[l]=Mean[n]; FWHM[l]=FWHM[n];
					Mean[n]=aux; FWHM[n]=aux2;}
			}
		}
		
		//linear regresion to get slope+offset (minimum squares)
		for(int l=0; l<4; l++){sumxy=0; sumxx=0; sumx=0; sumy=0;}
		for(int l=0; l<4; l++){		
			sumxy+=Mean[l]*energy[l];
			sumxx+=Mean[l]*Mean[l];
			sumx+=Mean[l]; sumy+=energy[l];
		}

		slope=(sumxy-sumx*sumy/4)/(sumxx-sumx*sumx/4);
		offset=(sumxx*sumy-sumx*sumxy)/(4*sumxx-sumx*sumx);
		//offset=(sumy-slope*sumx)/4;
		dynrange=4096*slope+offset;
	
		for(int l=0; l<4; l++) res[l]=FWHM[l]*slope;

		// write means and fwhms down in output file:
		//file << strip <<"	"<< Mean[0] <<"	"<< Mean[1] <<"	"<< Mean[2] <<"	"<< Mean[3] <<'\n';
		//file << strip <<"	"<< FWHM[0] <<"	"<< FWHM[1] <<"	"<< FWHM[2] <<"	"<< FWHM[3] <<'\n';

		// to see in terminal:
		/*cout << "offset: " << offset << " slope: " << slope << " res (keV @ 5804.96): " << res[3] 
										<< " Dyn. range (keV): " << dynrange <<'\n';*/

		// write the file
		file << offset << " " << slope << " " << ' ' << '\n';
		//for(int l=0; l<4; l++) 	file << res[l] << " " <<' ';
		//file << dynrange << '\n';		

	    }
	}

    //c1->Close();
    file.close();
}