BeEvent.h

/*
 * BeEvent.h
 *
 *  Created on: 9.12.2009
 *      Author: Vratislav
 */

#pragma once

#include <TObject.h>
#include <TFile.h>
#include <TTree.h>
#include <TString.h>
#include <TRandom3.h>
#include <TMath.h>
#include <TVector3.h>
#include <TRotation.h>
#include <TLorentzVector.h>
#include <TLorentzRotation.h>
#include <TCutG.h>

#include "../AculData/AculRaw.h"
#include "../AculData/AculCalibration.h"
#include "../Detectors/AnnularDetector.h"
#include "../TELoss/TELoss.h"
#include "../AculData/ConfigDictionary.h"

#include<iostream>

using 	std::cout;
using 	std::endl;

//TMath
using	TMath::Power;
using	TMath::Sqrt;

#define NOESEC 16	//number of sectors carrying full dE-E information
#define NORINGS 32	//number of all rings in the first detector of each telescope
#define NOHITS 16	//maximal number of hits
#define NOTEL 2		//number of telescopes

//materials used for simulations
//Si
#define SI_RHO 2.321 	//in g.cm^-3
#define SI_NELEMENTS 1  //number of elements in material
#define SI_1_A	28.086	//A of the first element
#define SI_1_Z	14.		//Z of the first element
#define SI_1_W	1.		//weight of the first element
//CsI
#define CSI_RHO 4.51 	//in g.cm^-3
#define CSI_NELEMENTS 2  //number of elements in material
#define CSI_1_A	132.91 	//A of the first element
#define CSI_1_Z	55.		//Z of the first element
#define CSI_1_W	0.5		//weight of the first element
#define CSI_2_A	126.9	//A of the second element
#define CSI_2_Z	53.		//Z of the second element
#define CSI_2_W	0.5		//weight of the second element

//target properties (19.8.2012)
#define T_NELEMENTS 1	//number of elements in material
#define T_1_A 1.008		//A of the first element
#define T_1_Z 1.		//Z of the first element
#define T_1_W 1.		//weight of the first element
#define T_NORMAL_RHO	0.00008895	//density in g/cm^3 at 273 K and 10^5 Pa
#define T_TEMPERATURE	30.			//temperature in K
#define T_PRESSURE		4.			//pressure in bars
#define T_THICKNESS		6000.		//target thickness in mcm

#define T_WIN_THICK		6.			//window thickness in mcm
#define T_WIN_RHO 8.
#define T_WIN_NELEMENTS 3 //number of elements in material
#define T_WIN_1_A 51.996		//A of the first element
#define T_WIN_1_Z 24.		//Z of the first element
#define T_WIN_1_W .08		//weight of the first element
#define T_WIN_2_A 55.847		//A of the first element
#define T_WIN_2_Z 26.		//Z of the first element
#define T_WIN_2_W .74		//weight of the first element
#define T_WIN_3_A 58.69		//A of the first element
#define T_WIN_3_Z 28.		//Z of the first element
#define T_WIN_3_W .18		//weight of the first element

//particles
#define kPROTON 11	//ZA
#define kALPHA 	24	//ZA

//masses
#define kPMASS	938.272		//in MeV/c^2
#define kNMASS	939.565		//in MeV/c^2
#define kAMASS	3728.400	//in MeV/c^2
#define	kLiMASS	5603.050	//in MeV/c^2
#define kBeMASS 5607.338	//in MeV/c^2

class BeEvent {

public:			//melo by byt private

	Int_t	fTrigger;		///also used for state identificator fTrigger == 1J for state J, i.e 10 for 0+, 12 for 0+, 11 for 1-
	Int_t	fn[NOESEC];

	Double_t	fCalCorr[NOESEC];

	//detectors
	AnnularDetector fT11;
	AnnularDetector fT12;
	AnnularDetector	fT13;
	AnnularDetector fT21;
	AnnularDetector fT22;
	AnnularDetector fT23;

	//deltaE in first detector, 32 sectors converted to 16 sectors
	Double_t	fdE1[NOESEC];
	Double_t	fdE2[NOESEC];

	//energy and multiplicity information from CsI detectors

	Double_t	fE13aSec[NOESEC];
	Double_t	fE13pSec[NOESEC];
	Double_t	fE23aSec[NOESEC];
	Double_t	fE23pSec[NOESEC];

	Int_t	fMult13a;
	Int_t	fMult23a;
	Int_t	fMult13p;
	Int_t	fMult23p;

	//time information from Tx2
	Double_t	fTau[NOTEL][NOESEC];	//time in ns

	//hits
	Int_t	fNOHits[NOTEL];	//general number of hits with dE-E, phi and theta information
	Int_t	fNOSHits[NOTEL];
	Int_t	fNORHits[NOTEL];	//hits in first detector of telescope detected only by rings
	Int_t	fNOaHits[NOTEL];	//alpha particle hits before choice by time
	Int_t	fNOpHits[NOTEL];	//proton hits before choice by time
	Int_t	fNOtHits[NOTEL];	//hits searched by time

	Int_t		fSN[NOTEL][NOHITS];
	Int_t		fTauSN[NOTEL][NOHITS];	//time hit sector number
	Int_t		fRN[NOTEL][NOHITS];
	Double_t	fE1Sec[NOTEL][NOHITS];	//hit energy in 1st layer (sectors) of particular telescope
	Double_t	fE1Ring[NOTEL][NOHITS];	//hit energy in 1st layer (rings) of particular telescope
	Double_t	fE2Sec[NOTEL][NOHITS];	//hit energy in 2nd layer of particular telescope
	Double_t	fE3Sec[NOTEL][NOHITS];	//hit energy in 3rd layer of particular telescope
	Double_t	fTheta[NOTEL][NOHITS];	//in rad for lab system
	Double_t	fPhi[NOTEL][NOHITS];	//in rad for lab system
	Double_t	fhTau[NOTEL][NOHITS];
	Double_t	fT[NOTEL][NOHITS];		//in MeV for lab system
	Double_t	fP[NOTEL][NOHITS];		//in MeV for ?? system
	Double_t	fdE1calc[NOTEL][NOHITS];	//calculated energy in the first detector, it is used for particle identification
//	Double_t	fE1[NOHITS];		//zatim bez pouziti
	Int_t		fID[NOTEL][NOHITS];			//particle identificator = ZA; proton==11, alpha==24
	Int_t		fPM[2*NOHITS];		//particle markers

//	//6Be information

	Int_t	fNOA;	//number of detected alphas related to triple coincidence Be decay
	Int_t	fNOP;	//number of detected protons related to triple coincidence Be decay
		//laboratory system
	//6Be
	TLorentzVector	f6BeLab;
	Double_t		f6BeThetaLab;	//Theta_lab in rad
	Double_t		f6BePhiLab;		//Phi_lab in rad
	Double_t	 	f6BePcLab;		//abs(pc) in MeV
	Double_t		f6BeIM;			//invariant mass in MeV
	//neutron
	TLorentzVector	fNLab;
	Double_t		fNThetaLab;		//Theta_lab in rad
	Double_t		fNPhiLab;		//Phi_lab in rad
	Double_t		fNPcLab;		//invariant mass in MeV
	//proton1
	TLorentzVector	fP1Lab;
	Double_t		fP1ThetaLab;	//Theta_lab in rad
	Double_t		fP1PhiLab;		//Phi_lab in rad
	Double_t		fP1PcLab;		//invariant mass in MeV
	//proton2
	TLorentzVector	fP2Lab;
	Double_t		fP2ThetaLab;	//Theta_lab in rad
	Double_t		fP2PhiLab;		//Phi_lab in rad
	Double_t		fP2PcLab;		//invariant mass in MeV
	//alpha
	TLorentzVector	fALab;
	Double_t		fAThetaLab;	//Theta_lab in rad
	Double_t		fAPhiLab;		//Phi_lab in rad
	Double_t		fAPcLab;		//invariant mass in MeV
	//protons relative kinetic energy
	Double_t		fPPTrel;

		//CM system 6Li-p
	//6Be
	TLorentzVector	f6BeCM1;
	Double_t		f6BeThetaCM1;	//Theta_cm in rad
	Double_t		f6BePhiCM1;		//Phi_lab in rad
	Double_t 		f6BePcCM1;		//abs(pc) in MeV

		//CM system 6Be
	//6Be
	TLorentzVector	f6BeCM;
	Double_t		f6BeThetaCM;	//Theta_cm in rad
	Double_t		f6BePhiCM;		//Phi_lab in rad
	Double_t 		f6BePcCM;		//abs(pc) in MeV
	//alpha
	TLorentzVector	fACM;
	Double_t		fAThetaCM;		//Theta_cm in rad
	Double_t		fAPhiCM;		//Phi_lab in rad
	Double_t 		fAPcCM;			//abs(pc) in MeV
	//proton1
	TLorentzVector	fP1CM;
	Double_t		fP1ThetaCM;		//Theta_cm in rad
	Double_t		fP1PhiCM;		//Phi_lab in rad
	Double_t 		fP1PcCM;		//abs(pc) in MeV
	//proton2
	TLorentzVector	fP2CM;
	Double_t		fP2ThetaCM;		//Theta_cm in rad
	Double_t		fP2PhiCM;		//Phi_lab in rad
	Double_t 		fP2PcCM;		//abs(pc) in MeV
	//neutron
	TLorentzVector	fNCM;
	Double_t		fNThetaCM;		//Theta_cm in rad
	Double_t		fNPhiCM;		//Phi_lab in rad
	Double_t 		fNPcCM;			//abs(pc) in MeV

		//general information
	Double_t	fQLiP;			//Q of the reaction 6Li+p-->6Be+n

		//Jacobi coordinates	//fixme check them
	//"T" system
	Double_t	fTpp;			//"T" system
	Double_t	fTapp;			//"T" system
	Double_t	fCosThetaTk;	//"T" system
	//"Y" system
	Double_t	fTap;			//"Y" system
	Double_t	fTpap;			//"Y" system
	Double_t	fCosThetaYk;	//"Y" system

	//variables used in Identification functions
	//theta angles
	Double_t	fTheta1;	//!
	Double_t	fQ1;		//		//todo apparently no information saved
	Double_t 	fT6Li2;				//todo apparently no information saved
	Double_t	fTheta2;	//!

	//computed energies
	Double_t	fE11aSec[NOESEC];		//computed energy in T11, fE11cSec = f(E12, x)
	Double_t	fE11pSec[NOESEC];		//compSetQLiP()uted energy in T11, fE11cSec = f(E12, x)
	Double_t	fE21aSec[NOESEC];		//computed energy in T21, fE21cSec = f(E22, x)
	Double_t	fE21pSec[NOESEC];		//computed energy in T21, fE21cSec = f(E22, x)
	Double_t	fE12aSec[NOESEC];		//computed energy in T12, fE12cSec = f(E13, x)
	Double_t	fE12pSec[NOESEC];		//computed energy in T12, fE12cSec = f(E13, x)
	Double_t	fE22aSec[NOESEC];		//computed energy in T22, fE21cSec = f(E23, x)
	Double_t	fE22pSec[NOESEC];		//computed energy in T22, fE21cSec = f(E23, x)


	//used in function DeltaEId. and DeltaECal. functions, used just for verification
	//these variables need not be be written into file
	//first telescope
//	Double_t	fX11;		//!	//draha v T11, pouze docasne
//	Double_t	fX12;		//!	//draha v T12, pouze docasne
	Double_t	fXd112;		//!	//draha v mrtve vrstve mezi T11 a T12
	Double_t	fXd123;		//!	//draha v mrtve vrstve mezi T12 a T13
	Int_t		fRN1;		//!
	Int_t		fSN1[NOHITS];		//!

	//second telescope
//	Double_t	fX21;		//!	//draha v T11, pouze docasne
//	Double_t	fX22;		//!	//draha v T12, pouze docasne
	Double_t	fXd212;		//!	//draha v mrtve vrstve mezi T11 a T12
	Double_t	fXd223;		//!	//draha v mrtve vrstve mezi T12 a T13
	Int_t		fRN2;		//!
	Int_t		fSN2;		//!

	//used for verification of deltaE calibration method for first telescope:
	Double_t	fE11v_0s[NOESEC];	//!
	Double_t	fE11v_0r;			//!
	Double_t	fE11vs[NOESEC];		//!
	Double_t	fE11vr;				//!
	Double_t	fE12v_0[NOESEC];	//!
	Double_t	fE12v[NOESEC];		//!

	//used for verification of deltaE calibration method for first telescope:
	Double_t	fE21v_0s[NOESEC];	//!
	Double_t	fE21v_0r;			//!
	Double_t	fE21vs[NOESEC];		//!
	Double_t	fE21vr;				//!
	Double_t	fE22v_0[NOESEC];	//!
	Double_t	fE22v[NOESEC];		//!

	//used for deltaE calibration of CsI
	//alhpas
	Double_t	fE12_0a[NOESEC];	//!
	Double_t	fE12a[NOESEC];		//!
	Double_t	fE13_0a[NOESEC];	//!
	//protons
	Double_t	fE12_0p[NOESEC];	//!
	Double_t	fE12p[NOESEC];		//!
	Double_t	fE13_0p[NOESEC];	//!
	//alhpas
	Double_t	fE22_0a[NOESEC];	//!
	Double_t	fE22a[NOESEC];		//!
	Double_t	fE23_0a[NOESEC];	//!
	//protons
	Double_t	fE22_0p[NOESEC];	//!
	Double_t	fE22p[NOESEC];		//!
	Double_t	fE23_0p[NOESEC];	//!





	static TRandom3	ranE;			//!
	static TRandom3 ranChoice;		//!
	TELoss		mSiAlphaV;		//!  //silicon
	TELoss		mSiAlpha;		//!  //silicon
	TELoss		mSiP;			//!  //silicon
	TELoss		mCsIAlpha;		//!  //CsI
	TELoss		mCsIP;			//!  //CsI
	TELoss		mH_P;			//!  //protons in hydrogen
	TELoss		mH_Alpha;		//!  //alpha in hydrogen
	TELoss		mWin_P;			//!  //protons in target window
	TELoss		mWin_Alpha;		//!  //alpha in target window

	//particle indentification cuts
	TFile	*fCuts;			//!
	TCutG	*cutAlpha1;		//!		1st telescope
	TCutG	*cutProton1;	//!
	TCutG	*cutAlpha2;		//!		2nd telescope
	TCutG	*cutProton2;	//!

	//calibration parameters
	AculCalibration *calSi;		//!		calibration parameters for Si detectors
	AculCalibration *calCsIa;	//!		alpha parameters for CsI detectors
	AculCalibration *calCsIp;	//!		proton parameters for CsI detectors
	TString sSiCal;				//!		Si calibration parameters
	TString sCsIprotonsCal;		//!		CsI calibration parameters for protons
	TString sCsIalphaCal;		//!		CsI calibration parameters for alphas

	//parameters
	TString fConfigFile;	//!

	//geometry used for data reconstruction
	Double_t fT1ExpPos;		//!	in mm
	Double_t fT2ExpPos;		//!	in mm

	//beam characteristics used in DA
	Double_t fTBeamRec;		//!	in AMeV		used for reconstruction //where is it? in the centre? on the target window?
//	Double_t fTBeamRes;		//!	in MeV, sigma for gaus
	Double_t fBeamX;		//!	in cm, x position of beam at target
	Double_t fBeamY;		//!	in cm, y position of beam at target
	Double_t fBeamX_sigma;		//!	in cm; x position sigma of beam at target
	Double_t fBeamY_sigma;		//!	in cm; x position sigma of beam at target

	//detector thicknesses in mcm
	Double_t fX11_FD;		//!
	Double_t fX11;			//!
	Double_t fX11_BD;		//!

	Double_t fX12_FD;		//!
	Double_t fX12;			//!
	Double_t fX12_BD;		//!

	//#define	XD13F	14.	//fixme original value		probably in Si equivalent
	Double_t fX13_FD;		//!

	Double_t fX21_FD;		//!
	Double_t fX21;			//!
	Double_t fX21_BD;		//!

	Double_t fX22_FD;		//!
	Double_t fX22;			//!
	Double_t fX22_BD;		//!
	//#define	XD23F	14.		//fixme original value		probably in Si equivalent
	Double_t fX23_FD;		//!

	Double_t fCorrT1_0_7_CsI_A;		//!	correction for calibration of alphas in T1 CsI detectors
	Double_t fCorrT1_8_15_CsI_A;	//!	correction for calibration of alphas in T1 CsI detectors
	Double_t fCorrT1_0_7_CsI_P;		//!	correction for calibration of protons in T1 CsI detectors
	Double_t fCorrT1_8_15_CsI_P;	//!	correction for calibration of protons in T1 CsI detectors

	//telescope 2
	Double_t fCorrT2_0_7_CsI_A;		//!	correction for calibration of alphas in T2 CsI detectors	//probably same as MS (17. 12. 2012)
	Double_t fCorrT2_8_15_CsI_A;	//!	correction for calibration of alphas in T2 CsI detectors
	Double_t fCorrT2_0_7_CsI_P;		//!	correction for calibration of protons in T2 CsI detectors
	Double_t fCorrT2_8_15_CsI_P;	//!	correction for calibration of protons in T2 CsI detectors

	//thresholds
	Double_t fRSTolerance;		//! in MeV, tolerance between energy in sectors and rings for one hit
	Double_t fProtonThr;		//! in MeV, upper threshold for protons

	TString fCutFile;			//! file containing cuts
	TString fCalibFilePath;		//!	path to files with calibration parameters

public:
	BeEvent();
	BeEvent(const char* configfile/*, const char* cutfile,*/ /*const TString calibfilepath*/);	//todo finish implementation of this constructor
								//change beam kinetic energy constant to a parameter
								//probably the data member of BeEvent and then
								//constructor will be in form of BeEvent(Double_t tbeam)
//	BeEvent(BeEvent &);
	virtual ~BeEvent();
	ClassDef(BeEvent, 1);

	void ReadConfigFile();

	void	FillEvent(AculRaw *rawevent);
	void	FillEventOld(AculRaw *rawevent, AculCalibration *calSi, AculCalibration *calA, AculCalibration *calP);

	//functions for work with particular detectors

	//opravit cislovani ringu v prvnim detektoru
	void	SetChannels(AculRaw *rawevent);	//fill all variables with information in channels
	void	SetSiDetEnergies(AculRaw *rawevent);	//fill all variables with energy information for T11, T12, T21, T22
	void	SetSiDetEnergiesOld(AculRaw *rawevent, AculCalibration *cal);	//fill all variables with energy information for T11, T12, T21, T22
	void	SetSiDetTimes(AculRaw *rawevent);
	void	SetSiDetTimesOld(AculRaw *rawevent, AculCalibration *cal);
	void	SetCsIalphaDetEnergies(AculRaw *rawevent);
	void	SetCsIalphaDetEnergiesOld(AculRaw *rawevent, AculCalibration *cal);
	void	SetCsIprotonDetEnergies(AculRaw *rawevent);	//zatim nepouzivat
	void	SetCsIprotonDetEnergiesOld(AculRaw *rawevent, AculCalibration *cal);	//zatim nepouzivat
	void	SetDeltaE();
	void	SetTheta1(Int_t ring);
	void	SetTheta2(Int_t ring);

	//functions for work with hits
	void	SetSecHitEnergiesSi(Int_t tel); 	//for first telescope notel==0, for second notel==1
	void	SetSecHitEnergiesSiOld(Int_t tel); 	//for first telescope notel==0, for second notel==1
	void	SetRingHitEnergies(Int_t tel); 	//for first telescope notel==0, for second notel==1
	void	SetHitsEnergies(Int_t tel);
	void 	SetEnergyHits(Int_t tel);
	void	SetThetas(Int_t tel);
	void	SetThetasOld(Int_t tel);
	void	SetThetasNew(Int_t tel);	//fixme redo this function, it should take into account beam position
	void	SetPhis(Int_t tel);
	void	BeParticleIdentificationNew(Int_t tel);
	void	BeParticleIdentification(Int_t tel);
	void	ParticleMarkers();
	void	SetTimeHits(Int_t tel);
	void	SetTimeHitsNew(Int_t tel);
	void	SetHitsOld();
	void	SetHits();
	void	SetT(Int_t tel);		//particle kinetic energy
	void	SetTNew(Int_t tel);		//particle kinetic energy
	Double_t CalcKinE(const Int_t telescope, const Int_t hit,
			const Double_t cosTheta, const Double_t fDeadCsI,
			const Double_t bDeadSi2, const Double_t fDeadSi2,
			const Double_t bDeadSi1, const Double_t fDeadSi1);
	void	SetP(Int_t tel);		//particle size of impuls

	//general functions
	void	InitDetectors();
	Bool_t	WriteCondition();
	void 	Reset();
	Bool_t	WriteConditionPure();

	//physics
	void	SetProductsLab();	//set kinematical information about 6Be and n in laboratory system
	void	SetProductsLabNew();	//set kinematical information about 6Be and n in laboratory system
	void	SetQLiP();
	void	SetProductsCM1();	//set kinematical information about 6Be and n in CM system 6Li-p
	void	SetProductsCM();	//set kinematical information about 6Be and n in CM system of 6Be

	//functions used for calibration
	void	Set6LiQ();
//	void	DeltaEIdentification1old();		//todo this function can be erased
//	void	DeltaEIdentification1();
//	void	DeltaEIdentification2();
	void	DeltaECalibration1();	//first telescope	//todo probably can be erased
	void	DeltaECalibration2();	//second telescope	//todo probably can be erased

private:
	Double_t	CsIEnergy(Int_t _address, Int_t _subaddress, AculRaw *_rawevent, AculCalibration *_cal, Double_t _x0, Double_t _ethr);
	Double_t	CsIcorrectedEnergy(Int_t _address, Int_t _subaddress, AculRaw *_rawevent, AculCalibration *_cal, Double_t _x0, Double_t _ethr);
	Double_t	CompareLosses(Double_t de_ex, Double_t e, Double_t x, Int_t particleID);  //de_ex experimental dE, e experimental E in second layer, x range in the dE layer, ID==24 for alpha, ID==11 for proton
	void		HitIdentification(const Int_t tel, const Int_t hit, const Double_t eCsI, const Int_t particle, const Double_t amin, const Double_t amax, const Double_t bmin, const Double_t bmax, const Double_t cmin, const Double_t cmax);
	Int_t		NOParticlesTauE(Int_t tel, Int_t particleID);	//number of particular particles which are detected as hit in proper time window
	Int_t		GetNOSecPhi(Int_t tel, Int_t esec);	//get number of physical sector in first telescope of each detector

	//CsI calibrations
	void		SetCalCorrs();
public:
	//relativistic kinematics
	static Double_t	PC2(Double_t T, Double_t mc2);			//neni potrebna, je obsazena v TLorentzVector
	static Double_t	T(Double_t pc2, Double_t mc2);			//v podstate taktez
	static Double_t	LawOfCosines(Double_t pb, Double_t pc, Double_t alpha);
	static Double_t	ReactionQ(Double_t Ta1, Double_t Ta2, Double_t mac2, Double_t mbc2, Double_t theta);
	static Double_t	ReactionQ(Double_t Ea, Double_t Eb, Double_t Ec, Double_t Ed);

	//auxiliary
	static Double_t	ReducedMass(Double_t m1, Double_t m2);
private:
	void ReadCuts(/*const char* cutfile*/);
	void CreateSiELosses();
	void CreateCsIELosses();
	void CreateTargetELosses();

};