AculCalibration.h

#pragma once

#include <TObject.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1I.h>
#include <TPolyMarker.h>
#include <TF1.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TGraph.h>
#include <TObjArray.h>
#include <TRandom3.h>
#include <THStack.h>
#include <TString.h>
#include <TSpectrum.h>

#include <iostream>
#include <fstream>
#include <iomanip>
#include <sstream>

#define DEFAULTNOPEAKS 20
#define	BLOCKSNUMBER 100
#define ADDRESSNUMBER 32

using std::cout;
using std::endl;
using std::setw;
using std::setprecision;
using std::stringstream;
using std::ostringstream;

class AculCalibration : public TObject
{


public:

	//smysl jako verejne globalni promenne maji:
	//fNOSpectra	- pocet zkalibrovanych spekter
	//????			- pocet spravne zkalibrovanych spekter
	//????			- pocet nespravne zkalibrovanych spekter
	//fEnergy[4]	- tabulka s energiemi piku, nacita se zvenci

//private:


//	TObjArray	*fHRawList;			//list of raw histograms, list is set to owner
	TObjArray	fHRawList;			//list of raw histograms, list is set to owner
	TObjArray	fHAnalyzedList;	//list of fitted and analyzed histograms, list is set to owner
	TObjArray	fHEnergyList;		//list of calibrated histograms, list is set to owner
	THStack		fHEnergyStack;		//some stack
	THStack		*fCurrentHStack;
	//dodelat current histograms
	TObjArray	fCurrentHistList;



//	TRandom3	*fRanGen;
//	static TRandom3 fRanGen;	//!
//	THStack		*fHStack;

	//parameters to be read from file
	Int_t		kRaNOPEAKS;
	Double_t	fEnergy[DEFAULTNOPEAKS];
	Double_t	fLowerChannel;
	Double_t	fUpperChannel;
	Double_t	fLowerPeakRelativeHight;	//pouziva se, private
	Double_t	fUpperPeakRelativeHight;	//pouziva se, private, nastavit nenulovou prednastavenou hodnotu
	Double_t	fPeakPositionTolerance;		//pouziva se, private
	Width_t		fFitFuncLineWidth;			//private
	Double_t	fFitMinSigma;				//pouziva se, private
	Double_t	fFitPeakThreshold;			//pouziva se, private, prozkoumat, k cemu vlastne slouzi ve fci ShowPeaks, popremyslet o vhodnem prednastaveni v konstruktoru

	//these variables are the main for the class
	Double_t	fAOld[BLOCKSNUMBER][ADDRESSNUMBER];	//calibration parameter, f(x) = fAOld*x + fBOld
	Double_t	fBOld[BLOCKSNUMBER][ADDRESSNUMBER];	//calibration parameter, f(x) = fAOld*x + fBOld
	Double_t	fC[BLOCKSNUMBER][ADDRESSNUMBER];	//treti kalibracni parametr, jine zavislosti nez pol1
	Double_t	fD[BLOCKSNUMBER][ADDRESSNUMBER];	//ctvrty kalibracni parametr

	TArrayD fA;
	TArrayD fB;




	//smysl je velmi pochybny
	TFile		*fCalInformation;

//private:
	Double_t	fPeak[DEFAULTNOPEAKS];				//v teto promenne je ulozena momentalni hodnota piku v kanalech, zejmena v jedne fci, mozno udelat ji jako lokalni, bude navratovou hodnotou fce PeaksFitting, predelat delku pole

public:
	AculCalibration();
	//default constructor
	AculCalibration(const char* parfile);		//
	virtual ~AculCalibration();

	Bool_t	SetInputParameters(const char* inputparfile = "parforcal.par");
	// Function which loads text file containing parameters for calibration
	//
	//	-inputparfile: file containing information on calibration source
	//	-In file with the data must be preserved systematic
	//	-There cannot be whitelines (probably)
	//
	//Example of "parforcal.par" file
	//................................................................
	//223Ra
	//
	//4		nopeaks		//number of peaks
	//4.415		E1		//in MeV
	//5.153		E2		//in MeV
	//5.683		E3		//in MeV
	//7.419		E4		//in MeV
	//100		lowerchannel	//in channels
	//4096		upperchannel	//in channels
	//0.1		lowerpeakhight		//in relative units
	//0.1		upperpeakhight		//in relative units
	//0.1		peakpositiontolerance	//in relative units
	//2		fitfunctionlinewidth	//integer 1 - 10
	//5		minfitsigma		//minimal sigma of the peaks to be fitted
	//0.4		fithightthreshold	//
	//................................................................

	Bool_t	SetCalibrationParameters(const char* calparfile);
	//Loads the file with calibration parameters
	//
	//  calparfile: file containing calibration parameters in format: crate number, address, subaddress, fAOld, fBOld, fC, fD
	//
	//	allowed comment characters: *, #, %, //

	void	PrintInputParameters();
	//I hope this is selfunderstanding function

	void	PrintCalibrationParameters(const Int_t blockmin = 1, const Int_t blockmax = BLOCKSNUMBER - 1);
	//Print calibration parameters fAOld, fBOld, fC, fD which are currently saved in object AculCalibration
	//
	//  blockmin: minimum  block data are displayed for
	//  blockmax: maximum  block data are displayed for

	Bool_t	CalculateCalibParameters(const char* inputfile, const char* block,
			const Int_t address, const char* treename, Int_t lowerchannel = 0,
			Int_t upperchannel = 4095, Int_t nEBins = 1000, Int_t lowersubaddress = 0,
			Int_t uppersubaddress = ADDRESSNUMBER-1); //calculate calibration parameters for given block in given file
	//function is not completely ready to use
	//
	//function for calculation of calibrate parameters for DAQ system based on "Go4"
	//
	//  inputfile: root file with calibration spectra
	//  block: block name to be calibrated
	//  lowerchannel: minimal channel from which the spectrum will be analysed
	//  upperchannel: maximal channel up to which the spectrum will be analysed
	//  lowersubaddress: block subaddress
	//  uppersubaddress: block subbaddress

	void CalibrateRawSpectra();

	UShort_t SomeFunction(UShort_t x, int len);

	template<typename T> T SomeFunction(T x, int len)
	{
		/*T max = x[0];
		for(int i = 1; i < len; i++)
			if(max < x[i])
				max = x[i];*/

//		return max;
		return x;
	}



	void CalibrateRawSpectra(const char* inputfile, const char* block, const Int_t address, const char* treename, Int_t lowerchannel = 0, Int_t upperchannel = 4095, Int_t nEBins = 1000, Int_t lowersubaddress = 0, Int_t uppersubaddress = ADDRESSNUMBER-1);

//	Bool_t	EnergyPositions(const char* inputfile, const char* block,
//			const Int_t address, const Int_t address, const char* treename, Int_t lowerchannel = 0,
//			Int_t upperchannel = 4095, Int_t lowersubaddress = 0, Int_t uppersubaddress = ADDRESSNUMBER-1);

	Int_t	PeaksFitting(TH1* hSpectrum, Option_t* option = "", Double_t sigmamin = 2);
	//
	//possible options: "V", "Q", ""

	Int_t	SearchPeaks(const TH1 *hin, Double_t sigma = 2, Option_t *option = "", const Int_t searchedpeaks = 100);

	void	FillRawSpectraFile(const char* rawdatafile, const char* block, const char* treename, TCanvas* rawCanvas = NULL, Option_t *option = "", Int_t xaxismin = 0, Int_t xaxismax = 4096);

	void	ShowRawSpectra(const char* filename, const Int_t block, TCanvas* rawCanvas = NULL, Int_t xaxismin = 0, Int_t xaxismax = 4096, /*TObjArray* histList = NULL,*/ const Int_t subaddress = 16);
	//this function is written for old format of raw data
	//do not use this function!!!!!!!!!!!
	void	ShowSpectra(const char* filename, TCanvas* rawCanvas = NULL, Option_t *option = "", Int_t xaxismin = 0, Int_t xaxismax = 4096, /*TObjArray* histList = NULL,*/ const Int_t subaddress = ADDRESSNUMBER);
	void	ShowAnalyzedSpectra(const char* filename, TCanvas* fittedRawCanvas = NULL, Int_t xaxismin = 0, Int_t xaxismax = 4096, Int_t subaddress = 16);
	//This function displays analyzed spectrum from a file, divides the canvas into a sufficient number of pads and displays
	//spectrums of each  block subadress on the suitable pads or displays one  selected spectrum .
	//Selects the peaks in the histogram and displays on the histogram how the spectrum were fitted.
	//
	//  filename: file .root containing analysed spectra
	//  fittedRawCanvas: canvas on which one you will see the spectrum
	//  xaxismin: Minimum channel, which will be displayed
	//  xaxismax: Maximum channel, which will be displayed
	//  subaddress:

	void	ShowEnergySpectra(const char *filename, TCanvas* energyCanvas = NULL, const Int_t subaddress = 16, Option_t* option = "", Double_t xaxismin = 0., Double_t xaxismax = 10.); //option: "sum", "c", "+",
	//Displays the spectrum  of the selected subbaddress block in MeV
	//
	//  filename: file .root containing calibrated spectra in MeV
	//  energyCanvas: : canvas on which one you will see the spectrum
	//  subaddress: block subaddress  which will be drawn
	//  option: sum ,+ ,c
	//  xaxismin: Minimum channel, which will be displayed
	//  xaxismax: Maximum channel, which will be displayed

	void DivideCanvas(TCanvas *c1, Int_t inputs);



	//dodelat funkce TTree* Get...(...) pro raw, anal i E spectra
	Bool_t	AddCalFileToList(const char* calfilelist = "CalFileList.log");

	void	ClearHistograms(Option_t* option = "");
	void	DeleteStacks(Option_t* option = "");

	void	MakeCalibrationFile(Char_t* calibrationfile, Char_t* calfilelist);

	void	Reset();

//private:

	ClassDef(AculCalibration, 1);

};