fillChain.cxx

#if !defined(__CLING__)

#include "TSystem.h"
#include "TFile.h"
#include "TTree.h"
#include "TVector3.h"
#include "TH2F.h"
#include "../src/TNeEvent.h"
#include "../src/TNeDet16.h"
#include "TStopwatch.h"

#endif

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

Int_t GetClustersMWPC(unsigned short n, unsigned short *x)
{
	Int_t clusters = 1;

	if(n == 0) clusters = 0.;
	else if(n > 1 && n<=32)
	{
		for(Int_t i = 1; i < n; ++i)
		{
//			if (x[i]>50 || x[i-1]>50) continue;
			if( (abs(x[i] - x[i-1]) > 1) && (abs(x[i] - x[i-1]) < 32) ) clusters += 1.;
		}
	}
	return clusters;
} //--------------------------------------------------------------------

Double_t GetClusterPositionMWPC(unsigned short n, unsigned short *x, Float_t wireStep, Float_t planeOffset=0.)
{
	Double_t position = -100.;




//	if(n == 0) clusters = 0.;
	if(n > 0 && n<=32)
	{

//		x2p = (revent->x2[0] + 0.5 - 16)*kMWPCwireStepX2 + MWPC2_X_offset;
		position = (x[0]+n/2. + 0.5 -16)*wireStep + planeOffset;
//		cout << n << endl;
	}
	return position;
} //--------------------------------------------------------------------

void fillTree(const TString beam = "he", Int_t nofile = 0, const Int_t noevents = 0) {

	TString inFile;
	if (beam.Contains("he")) inFile.Form("~/data/exp1804/h5_14_%02d.root", nofile);
	if (beam.Contains("be")) inFile.Form("~/data/exp1804/be10_03_%d0.root", nofile);		//files 00,10,...,90
																//where 70 is run 01
																//		80 is run 02
																//		90 is run 05
//	inFile.Form("~/data/exp1804/calib/si_20_03.root");

	TString outFile;
	if (beam.Contains("he")) outFile.Form("~/data/exp1804/h5_14_%02d_calib.root", nofile);
	if (beam.Contains("be")) outFile.Form("~/data/exp1804/be10_03_%d0_calib.root", nofile);	//files 00,10,...,60
																	//where 70 is run 01
																	//		80 is run 02
																	//		90 is run 05
//	outFile.Form("~/data/exp1804/calib/si_20_03_calib.root");

	cout << "Input file: " << inFile << endl;
	cout << "Output file: " << outFile << endl;

	///////////////////////////////////////////////////
	// Input file initialization
	///////////////////////////////////////////////////

	TFile *fr = new TFile(inFile);
	TTree *tr = (TTree*)fr->Get("AnalysisxTree");

	TNeEvent *revent = new TNeEvent();
	tr->SetBranchAddress("NeEvent.", &revent);


	///////////////////////////////////////////////////
	// Output file initialization
	///////////////////////////////////////////////////

	TFile *fw = new TFile(outFile, "RECREATE");
	TTree *tw = new TTree("cal", "Calibrated information");

	///////////////////////////////////////////////////
	// Output tree variables
	///////////////////////////////////////////////////

	Int_t trigger;

	Float_t SQ20E[16];
	Float_t SQ20Ecorr[16];
	Float_t SQ20EcorrHit[16];
	Float_t SQ20EcorrHitC[16];
	Float_t SQ20Esum;

	Float_t SQ20time[16];
	Int_t SQ20timeMult;
	Float_t SQLXtime[32];
	Int_t SQLXtimeMult;
//	Float_t SQLXtimeSum;

	Float_t SQLXE[32];
	Float_t SQLXEsum;

	Float_t SQLXEtimeFiltered[32];
	Float_t SQLXEtimeFilteredSum;

	Float_t SQLYE[16];
	Float_t SQLYEsum;
	Float_t SQLYtime[16];
	Int_t SQLYtimeMult;

	Float_t SQLYEtimeFiltered[16];
	Float_t SQLYEtimeFilteredSum;

	Int_t SQLXmult;
	Int_t SQLYmult;
	Int_t SQ20EcorrMult;
	Int_t SQ20EcorrHitMult;
	Int_t SQ20EcorrHitCMult;

	Bool_t CsI_L_veto;
	Int_t expectedThinStrip;


	Float_t SQRXE[32];
	Float_t SQRXEsum;
	Int_t SQRXmult;
	Float_t TOF, dEbeam;

	Float_t tF5;

	//MWPC, wire multiplicity
	Float_t x1p, y1p, x2p, y2p, xt, yt;

	//MWPC, cluster multiplicity
	Int_t x1MultC, x2MultC, y1MultC, y2MultC;
	Float_t x1c, y1c, x2c, y2c, xtc, ytc;


	//left 1 mm position
	Float_t x1mm, y1mm;

	Float_t xThin, yThin;

	Int_t mapXbin, mapYbin;

	TVector3 vHit1mm;
	TVector3 vNorm(0.,0.,1.);
	Double_t angleLeft;

	///////////////////////////////////////////////////
	// Output tree branches initialization
	///////////////////////////////////////////////////

	tw->Branch("trigger",&trigger,"trigger/I");

	tw->Branch("angleLeft",&angleLeft,"angleLeft/D");
	tw->Branch("x1mm",&x1mm,"x1mm/F");
	tw->Branch("y1mm",&y1mm,"y1mm/F");
	tw->Branch("xThin",&xThin,"xThin/F");
	tw->Branch("yThin",&yThin,"yThin/F");

	tw->Branch("mapXbin",&mapXbin,"mapXbin/I");
	tw->Branch("mapYbin",&mapYbin,"mapYbin/I");

	tw->Branch("SQ20E",SQ20E,"SQ20E[16]/F");
	tw->Branch("SQ20Ecorr",SQ20Ecorr,"SQ20Ecorr[16]/F");
	tw->Branch("SQ20EcorrHit",SQ20EcorrHit,"SQ20EcorrHit[16]/F");
	tw->Branch("SQ20EcorrHitC",SQ20EcorrHitC,"SQ20EcorrHitC[16]/F");
	tw->Branch("SQ20Esum",&SQ20Esum,"SQ20Esum/F");

	tw->Branch("SQ20time",SQ20time,"SQ20time[16]/F");
	tw->Branch("SQ20timeMult",&SQ20timeMult,"SQ20timeMult/I");
//	tw->Branch("SQLXtimeSum",&SQLXtimeSum,"SQLXtimeSum/F");

	tw->Branch("SQLXtime",SQLXtime,"SQLXtime[32]/F");
	tw->Branch("SQLXtimeMult",&SQLXtimeMult,"SQLXtimeMult/I");

	tw->Branch("SQLYtime",SQLYtime,"SQLYtime[16]/F");
	tw->Branch("SQLYtimeMult",&SQLYtimeMult,"SQLYtimeMult/I");

	tw->Branch("SQLXE",SQLXE,"SQLXE[32]/F");
	tw->Branch("SQLXEsum",&SQLXEsum,"SQLXE/F");

	tw->Branch("SQLXEtimeFiltered",SQLXEtimeFiltered,"SQLXEtimeFiltered[32]/F");
	tw->Branch("SQLXEsumtimeFilteredSum",&SQLXEtimeFilteredSum,"SQLXEtimeFilteredSum/F");

	tw->Branch("SQLYEtimeFiltered",SQLYEtimeFiltered,"SQLYEtimeFiltered[16]/F");
	tw->Branch("SQLYEsumtimeFilteredSum",&SQLYEtimeFilteredSum,"SQLYEtimeFilteredSum/F");

	tw->Branch("SQLYE",SQLYE,"SQLYE[16]/F");
	tw->Branch("SQLYEsum",&SQLYEsum,"SQLYEsum/F");

	tw->Branch("SQLXmult",&SQLXmult,"SQLXmult/I");
	tw->Branch("SQLYmult",&SQLYmult,"SQLYmult/I");
	tw->Branch("SQ20EcorrMult",&SQ20EcorrMult,"SQ20EcorrMult/I");
	tw->Branch("SQ20EcorrHitCMult",&SQ20EcorrHitCMult,"SQ20EcorrHitCMult/I");
//	tw->Branch("SQ20EcorrMult",&SQ20EcorrMult,"SQ20EcorrMult/I");

	tw->Branch("CsI_L_veto", &CsI_L_veto, "CsI_L_veto/O");
	tw->Branch("expectedThinStrip",&expectedThinStrip,"expectedThinStrip/I");

	tw->Branch("SQRXE",SQRXE,"SQRXE[32]/F");
	tw->Branch("SQRXEsum",&SQRXEsum,"SQRXE/F");
	tw->Branch("SQRXmult",&SQRXmult,"SQRXmult/I");

	tw->Branch("TOF", &TOF, "TOF/F");
	tw->Branch("dEbeam", &dEbeam, "dEbeam/F");
	tw->Branch("tF5", &tF5, "tF5/F");

	tw->Branch("x1p", &x1p, "x1p/F");
	tw->Branch("y1p", &y1p, "y1p/F");
	tw->Branch("x2p", &x2p, "x2p/F");
	tw->Branch("y2p", &y2p, "y2p/F");

	tw->Branch("xt", &xt, "xt/F");
	tw->Branch("yt", &yt, "yt/F");

	tw->Branch("x1MultC",&x1MultC,"x1MultC/I");
	tw->Branch("x2MultC",&x2MultC,"x2MultC/I");
	tw->Branch("y1MultC",&y1MultC,"y1MultC/I");
	tw->Branch("y2MultC",&y2MultC,"y2MultC/I");

	tw->Branch("x1c", &x1c, "x1c/F");
	tw->Branch("y1c", &y1c, "y1c/F");
	tw->Branch("x2c", &x2c, "x2c/F");
	tw->Branch("y2c", &y2c, "y2c/F");

	tw->Branch("xtc", &xtc, "xtc/F");
	tw->Branch("ytc", &ytc, "ytc/F");


	Long64_t nevents;
	if (noevents == 0) nevents = tr->GetEntries();
	else nevents = noevents;
	if (nevents > tr->GetEntries()) nevents = tr->GetEntries();

	///////////////////////////////////////////////////
	// calibration coefficients and thickness map
	///////////////////////////////////////////////////


//	TNeDet16 *pSQX_L_EC = new TNeDet16("SQX_L_EC");
//	TNeDet16 pSQX_L_EC("../SQX_L_EC");
	TNeDet16 pSQX_L_EC("./SQX_L");
	pSQX_L_EC.ReadData();

	TNeDet16 pSQX_L_TC("./SQX_L_time");
	pSQX_L_TC.ReadData();

	TNeDet16 pSQY_L_EC("./SQY_L");
	pSQY_L_EC.ReadData();

	TNeDet16 pSQX_R_EC("../SQX_R_EC");
	pSQX_R_EC.ReadData();

	TNeDet16 pSQ20_EC("./sq20_58");
	pSQ20_EC.ReadData();

	Float_t energy = 0;

	cout << nevents << " entries will be treated." << endl;

//	TFile fThickness("thicknessPreliminary.root", "OPEN");
	TFile fThickness("thicknessFinal.root", "OPEN");
	fr->cd();
	TH2F *hThickness = new TH2F(*(TH2F*)fThickness.Get("hTh"));
//	hThickness->Draw("col");


	///////////////////////////////////////////////////
	// Parameters related to geometry
	///////////////////////////////////////////////////

	//all distances in mm

	//MWPC
	const Float_t kMWPCz1 = -815.;	//z coordinate of the center of MWPC1
	const Float_t kMWPCz2 = -270.;	//z coordinate of the center of MWPC2

	//step between two wires;
	//sign mean onward (+) or rearward (-) to physical axis
	const Float_t kMWPCwireStepX1 = -1.25;
	const Float_t kMWPCwireStepY1 = 1.25;		//step between two wires
	const Float_t kMWPCwireStepX2 = -1.25;		//step between two wires
	const Float_t kMWPCwireStepY2 = 1.25;		//step between two wires

	//offsets taken from S. Krupko
	const Float_t MWPC1_X_offset = -1.19;
	const Float_t MWPC1_Y_offset = -2.12;
	const Float_t MWPC2_X_offset = 0.2;
	const Float_t MWPC2_Y_offset = -1.02;

	//left telescope
	const Float_t z1mm = 230.;
	const Float_t zThin = 230.-53.6;
	//TODO: check signs and values
//	const Float_t xThinOffset = -3.;
//	const Float_t yThinOffset = -1.8;


	//TODO: add sign taking into account direction of numbering
	const Int_t kSQL_X_strips = -32;
	const Int_t kSQL_Y_strips = 16;
	const Int_t kSQL_20_strips = 16;

	//1mm left size
	const Float_t kSQL_X_size = 60.;
	const Float_t kSQL_Y_size = 60.;

	//thin detector
	const Double_t kSQ20_norm_thickness = 20.;
	const Float_t kSQ20_X_size = 50.;
	const Float_t kSQ20_Y_size = 50.;

	const Double_t kThinXoffset = -1.;
	const Double_t kThinYoffset = -1.8;


	//left CsI detectors
	const UShort_t CsIleftThr = 180;

	///////////////////////////////////////////////////
	// Individual thresholds
	///////////////////////////////////////////////////


	Float_t timeCorr[16] = {0, -2, 0, 2,
							-2.5, 5, -6, 4,
							7, 3, 3, 4,
							6, 4, -4, 13
	};

	Int_t thinChThresholds[16] = {
		113, 121, 125, 125,
		125, 132, 130, 115,
		118, 128, 132, 130,
		130, 130, 130, 112
	};

	Int_t thickLYChThresholds[16] = {
		4000, 130, 82, 110,
		87, 112, 82, 108,
		82, 112, 82, 108,
		80, 110, 90, 105,
	};

	Int_t thickLXChThresholds[32] = {
		60, 45, 45, 45,
		45, 4000, 4000, 45,
		45, 45, 45, 45,
		45, 45, 45, 45,
		/////
		45, 45, 45, 45,
		45, 45, 45, 45,
		45, 45, 45, 45,
		45, 45, 45, 60
	};

	Int_t minTimeLthin[16] {
		410, 410, 420, 415,
		420, 410, 415, 415,
		420, 420, 420, 420,
		410, 425, 0, 0
	};

	Int_t maxTimeLthin[16] {
		520, 525, 520, 530,
		520, 530, 525, 525,
		530, 530, 530, 520,
		510, 525, 0, 0
	};

	Double_t minTimeLX[32] {
		0, 0, 325, 328,
		330, 0, 0, 330,
		326, 328, 326, 335,
		326, 328, 328, 330,
		//
		420, 425, 420, 420,
		420, 420, 420, 420,
		430, 425, 420, 420,
		415, 415, 425, 415
	};

	Double_t maxTimeLX[32] {
		0, 0, 345, 345,
		350, 0, 0, 350,
		345, 345, 345, 350,
		340, 342, 345, 345,
		//
		460, 460, 460, 460,
		460, 460, 460, 460,
		460, 460, 460, 460,
		460, 460, 460, 453
	};

	fw->cd();

	//////////////////////////////////
	//event processing
	//////////////////////////////////

	for (Int_t eventNo = 0; eventNo < nevents; eventNo++) {

		tr->GetEvent(eventNo);

		trigger = revent->trigger;

		if (trigger != 3) continue;

		//reset values to be written

		SQ20Esum = 0.;
		SQLXEsum = 0.;
		SQLYEsum = 0.;
		SQRXEsum = 0.;
		SQRXmult = 0;

		SQLXEtimeFilteredSum = 0.;
		SQLYEtimeFilteredSum = 0.;

		SQLXmult = 0;
		SQLYmult = 0;
		SQ20EcorrMult = 0;
		SQ20EcorrHitMult = 0;
		SQ20EcorrHitCMult = 0;

		SQ20timeMult = 0;
		SQLXtimeMult = 0;
		SQLYtimeMult = 0;

		CsI_L_veto = kFALSE;
		expectedThinStrip = -1;

		x1p = -100.;
		y1p = -100.;
		x2p = -100.;
		y2p = -100.;
		xt = -100.;
		yt = -100.;

		x1c = -100.;
		y1c = -100.;
		x2c = -100.;
		y2c = -100.;
		xtc = -100.;
		ytc = -100.;

		x1MultC = 0;
		y1MultC = 0;
		x2MultC = 0;
		y2MultC = 0;


		///////////////////////////////////////////////
		//beam diagnostics
		///////////////////////////////////////////////

			///////////////////////////////////////////////
			//TOF
			///////////////////////////////////////////////

		dEbeam = 0.;
		TOF = 0.;
		if (revent->tF5[0]>0 && revent->tF5[1]>0 && revent->tF5[2]>0 && revent->tF5[3]>0
				&& revent->tF3[0]>0 && revent->tF3[1]>0 && revent->tF3[2]>0 && revent->tF3[3]>0
				&& revent->F5[0]>0 && revent->F5[1]>0 && revent->F5[2]>0 && revent->F5[3]>0
//				&& ( (revent->tF5[0]+revent->tF5[1]+revent->tF5[2]+revent->tF5[3]) - (revent->tF3[0]+revent->tF3[1]+revent->tF3[2]+revent->tF3[3]) )/4.*0.125+89.165<200
//				&& ( (tF5[0]+tF5[1]+tF5[2]+tF5[3]) - (tF3[0]+tF3[1]+tF3[2]+tF3[3]) )/4.*0.125+89.165>100
//				&& (F5[0]+F5[1]+F5[2]+F5[3])/4. < 2500
		) {
			dEbeam = (revent->F5[0]+revent->F5[1]+revent->F5[2]+revent->F5[3])/4.;
			TOF = ( (revent->tF5[0]+revent->tF5[1]+revent->tF5[2]+revent->tF5[3]) - (revent->tF3[0]+revent->tF3[1]+revent->tF3[2]+revent->tF3[3]) )/4.*0.125+89.165;
		}

		if (TOF<166. || TOF>181.) continue;

			///////////////////////////////////////////////
			//MWPC's
			///////////////////////////////////////////////


		Bool_t flagMWPC = 1;

		tF5 = 0.125*(revent->tF5[0]+revent->tF5[1]+revent->tF5[2]+revent->tF5[3])/4.;
//		tF5 = 0.125*revent->tF5[0];

		//check for MWPC timing
		//TODO: parametrization of time thresholds
		if(revent->nx1>10 || revent->ny1>10 || revent->nx2>10 || revent->ny2>10) flagMWPC=0;
		if((0.125*revent->tMWPC[0]-tF5)<60 || (0.125*revent->tMWPC[0]-tF5)>77) flagMWPC=0;
		if((0.125*revent->tMWPC[1]-tF5)<60 || (0.125*revent->tMWPC[1]-tF5)>80) flagMWPC=0;
		if((0.125*revent->tMWPC[2]-tF5)<70 || (0.125*revent->tMWPC[2]-tF5)>90) flagMWPC=0;
		if((0.125*revent->tMWPC[3]-tF5)<60 || (0.125*revent->tMWPC[3]-tF5)>80) flagMWPC=0;

		if (flagMWPC==0) continue;

		//MWPC: work with wire multiplicity equal to 1
		if (flagMWPC!=0) {
			if (revent->x1[0]<1000 && revent->y1[0]<1000 && revent->nx1==1 && revent->ny1==1) {
//				x1p = (revent->x1[0]+0.5)*1.25-20. + MWPC1_X_offset;
//				y1p = (revent->y1[0]+0.5)*1.25-20. + MWPC1_Y_offset;

				x1p = (revent->x1[0] + 0.5 - 16)*kMWPCwireStepX1 + MWPC1_X_offset;
				y1p = (revent->y1[0] + 0.5 - 16)*kMWPCwireStepY1 + MWPC1_Y_offset;
			}

			if (revent->x2[0]<1000 && revent->y2[0]<1000 && revent->nx2==1 && revent->ny2==1) {
//				x2p = (revent->x2[0]+0.5)*1.25-20. + MWPC2_X_offset;
//				y2p = (revent->y2[0]+0.5)*1.25-20. + MWPC2_Y_offset;

				x2p = (revent->x2[0] + 0.5 - 16)*kMWPCwireStepX2 + MWPC2_X_offset;
				y2p = (revent->y2[0] + 0.5 - 16)*kMWPCwireStepY2 + MWPC2_Y_offset;
			}

			if (x1p > -50. && y1p > -50. && x2p > -50. && y2p > -50.) {
//				xt = x1p - (x1p - x2p)*((l12+lt)/l12);
//				yt = y1p - (y1p - y2p)*((l12+lt)/l12);
//				xt = (l12*x1p + (l12+lt)*(x2p-x1p))/(l12 - (x2p-x1p)*TMath::Tan(TMath::DegToRad()*12.));
//				yt = (y2p-y1p)*(xt-x1p)/(x2p-x1p) + y1p;

				xt = x1p - (x2p -x1p)*kMWPCz1/(kMWPCz2-kMWPCz1);
				yt = y1p - (y2p -y1p)*kMWPCz1/(kMWPCz2-kMWPCz1);

			}
		}//if wire multiplicity == 1

//		cout << flagMWPC << endl;

		//MWPC: work with cluster multiplicity equal to 1
		if (flagMWPC!=0) {
			x1MultC = GetClustersMWPC(revent->nx1, revent->x1);
			x2MultC = GetClustersMWPC(revent->nx2, revent->x2);
			y1MultC = GetClustersMWPC(revent->ny1, revent->y1);
			y2MultC = GetClustersMWPC(revent->ny2, revent->y2);

			//set wire step as input parameter for the function
			if (x1MultC==1 && y1MultC==1) {
				x1c = GetClusterPositionMWPC(revent->nx1, revent->x1, kMWPCwireStepX1, MWPC1_X_offset);
				y1c = GetClusterPositionMWPC(revent->ny1, revent->y1, kMWPCwireStepY1, MWPC1_Y_offset);
			}

			if (x2MultC==1 && y2MultC==1) {
				x2c = GetClusterPositionMWPC(revent->nx2, revent->x2, kMWPCwireStepX2, MWPC2_X_offset);
				y2c = GetClusterPositionMWPC(revent->ny2, revent->y2, kMWPCwireStepY2, MWPC2_Y_offset);
			}

			if (x1c > -50. && y1c > -50. && x2c > -50. && y2c > -50.) {
//				xtc = x1c - (x1c - x2c)*((l12+lt)/l12);
//				ytc = y1c - (y1c - y2c)*((l12+lt)/l12);
//				xtc = (l12*x1c + (l12+lt)*(x2c-x1c))/(l12 - (x2c-x1c)*TMath::Tan(TMath::DegToRad()*12.));
//				ytc = (y2c-y1c)*(xtc-x1c)/(x2c-x1c) + y1c;

				xt = x1c - (x2c -x1c)*kMWPCz1/(kMWPCz2-kMWPCz1);
				yt = y1c - (y2c -y1c)*kMWPCz1/(kMWPCz2-kMWPCz1);
			}
		}


		if (x1MultC!=1 || y1MultC!=1 || x2MultC!=1 || y2MultC!=1) continue;



		///////////////////////////////////////////////
		//work with Si detectors
		///////////////////////////////////////////////

//cout << endl;

		//CdI left -veto:
		for(Int_t i = 0; i<16; i++) {
			if (revent->CsI_L[i]>CsIleftThr) CsI_L_veto = kTRUE;
		}


		for (Int_t i = 0; i < 32; i++) {
//			cout << pSQX_L_EC.Energy(1, i) << endl;
			SQLXE[i] = 0;
			SQRXE[i] = 0;
			SQLXtime[i] = 0;
			SQLXEtimeFiltered[i] = 0.;

			if (revent->SQX_L[i] > thickLXChThresholds[i]) {
				energy = pSQX_L_EC.Energy(revent->SQX_L[i], i);
				SQLXE[i] = energy;
				SQLXEsum += SQLXE[i];
				SQLXmult++;
			}

			if (i<16 && i!=0 && i!=1
					&& SQLXE[i]>0
					&& revent->tSQX_L[i]*0.3>minTimeLX[i] && revent->tSQX_L[i]*0.3<maxTimeLX[i]) {
				SQLXtime[i] = revent->tSQX_L[i]*0.3;
				SQLXtimeMult++;
//				SQLXtimeSum = SQLXtimeSum + SQLXtime[i];
				SQLXEtimeFiltered[i]=energy;
				SQLXEtimeFilteredSum = SQLXEtimeFilteredSum + SQLXEtimeFiltered[i];
//				cout << SQLXtime[i] << endl;
//				cout << i << endl;
			}
			if (i>15
					&& SQLXE[i]>0
					&& revent->tSQX_L[i]*0.3>minTimeLX[i] && revent->tSQX_L[i]*0.3<maxTimeLX[i]) {
				SQLXtime[i] = revent->tSQX_L[i]*0.3;
				SQLXtimeMult++;
//				SQLXtimeSum = SQLXtimeSum + SQLXtime[i];
				SQLXEtimeFiltered[i] = energy;
				SQLXEtimeFilteredSum = SQLXEtimeFilteredSum + SQLXEtimeFiltered[i];
////				cout << SQLXtime[i] << endl;
			}

			energy = pSQX_R_EC.Energy(revent->SQX_R[i], i);

			if (energy > 1.0) {
				SQRXE[i] = energy;
				SQRXEsum += SQRXE[i];
				SQRXmult++;
			}
		}// for 32

		for (Int_t i = 0; i < 16; i++) {
			//			cout << pSQX_L_EC.Energy(1, i) << endl;

			SQLYE[i] = 0.;
			SQ20E[i] = 0.;
			SQ20time[i] = 0.;
			SQ20Ecorr[i] = 0.;
			SQ20EcorrHit[i] = 0.;
			SQ20EcorrHitC[i] = 0.;
			//				SQRXE[i] = 0;

			SQLYtime[i] = 0.;
			SQLYEtimeFiltered[i] = 0.;

			if (revent->SQ20[i] > thinChThresholds[i]) {
				energy = pSQ20_EC.Energy(revent->SQ20[i], i);
				SQ20E[i] = energy;
				SQ20Esum += SQ20E[i];
//				cout << energy << endl;
//				cout << i << "\t" << energy << "\t" << pSQ20_EC.a[i][0] << "\t" << pSQ20_EC.a[i][1] << "\t" << pSQ20_EC.a[i][2] << endl;
			}

			SQ20time[i] = revent->tSQ20[i]*0.3;
			if (SQ20E[i]>0 && SQ20time[i]>minTimeLthin[i] && SQ20time[i]<maxTimeLthin[i]) SQ20timeMult++;


//			if (i==0) continue;



			//				cout << energy << "\t" << pSQY_L_EC.a[i][0] << "\t" << pSQY_L_EC.a[i][1] << "\t" << pSQY_L_EC.a[i][2] << endl;
			//			cout << energy << endl;

			if (revent->SQY_L[i] > thickLYChThresholds[i]) {
				energy = pSQY_L_EC.Energy(revent->SQY_L[i], i);
				SQLYE[i] = energy;
				SQLYEsum += SQLYE[i];
				SQLYmult++;
			}

			Double_t time = revent->tSQY_L[i]*0.3 + timeCorr[i];

			if (SQLYE[i]>0 && time>325 && time<333) {
				SQLYtime[i] = time;
				SQLYtimeMult++;
				SQLYEtimeFiltered[i]=energy;
				SQLYEtimeFilteredSum = SQLYEtimeFilteredSum + SQLYEtimeFiltered[i];
			}

		}//for 16

//		cout << SQLXtimeMult << endl;

//		cout << endl;

		///////////////////////////////////////////////
		// correction for thickness inhomogeneity
		// when not taking into account hitting
		// particle tracking
		//
		// note: correction for thickness is not affected
		// by used system of coordinates
		///////////////////////////////////////////////

		Double_t currThickness;

		if (SQLXmult == 1 && SQLYmult == 1) {

			for (Int_t yi = 0; yi < kSQL_Y_strips; yi++) {

				//if energy in 1mm Y strip is reasonable
				if (SQLYE[yi]>0) {

					for (Int_t xi = 0; xi < TMath::Abs(kSQL_X_strips); xi++) {

						//if energy in 1mm X strip is reasonable
						if (SQLXE[xi] > 0) {
							currThickness = hThickness->GetBinContent(xi+1, yi+1);

							for (Int_t x20 = 0; x20 < kSQL_20_strips; x20++) {
								if (SQ20E[x20] > 0/*SQ20timeMult==1*/ && currThickness<30.) {

									SQ20Ecorr[x20] = (kSQ20_norm_thickness/currThickness)*SQ20E[x20];
									SQ20EcorrMult++;

								}
							}//for thin detector X strips
						}//if 1 mm X

					}//for thick detector X strips
				}//if 1 mm Y
			}//for thick detector Y strips

		}//if for correction according energy multiplicity




		///////////////////////////////////////////////
		// position of hit in left telescope
		///////////////////////////////////////////////

		//reset
		mapXbin = -1;
		mapYbin = -1;
		x1mm = -100.;
		y1mm = -100.;
		xThin = -100.;
		yThin = -100.;
		vHit1mm.SetXYZ(x1mm, y1mm, -1.);
		angleLeft = TMath::Pi();

		if (SQLXtimeMult==1) {
			for (Int_t i = 0; i<32; i++) {
				if (SQLXEtimeFiltered[i]>0) {
//					x1mm = -30. + (i+1./2.)*60./32.;
					x1mm = (i+0.5)*kSQL_X_size/kSQL_X_strips + kSQL_X_size/2.;	//ok
				}
			}
		}
//		cout << SQLYtimeMult << endl;

		if (SQLYtimeMult==1) {
			for (Int_t i = 0; i<16; i++) {
				if (SQLYEtimeFiltered[i]>0) {
//					y1mm = -30. + (i+1./2.)*60./16.;
					y1mm = -kSQL_Y_size/2. + (i+0.5)*kSQL_Y_size/kSQL_Y_strips;	//ok
				}
			}
		}

		//for MWPC wire multiplicity == 1

		//vHit1mm conversion to plane of thin detector
		if (SQLXtimeMult==1 && SQLYtimeMult==1 && xt>-50. && yt>-50.) {
			vHit1mm.SetXYZ(x1mm-xt, y1mm-yt, z1mm);
			vHit1mm *= zThin/z1mm;
			angleLeft = vHit1mm.Angle(vNorm);
		}


		if (SQLXtimeMult==1 && SQLYtimeMult==1
				//check if the particle passed through sensitive area of thin detector
				&& vHit1mm.X()+xt < (kSQ20_X_size/2. + kThinXoffset)
				&& vHit1mm.X()+xt > (-kSQ20_X_size/2. + kThinXoffset)
				&& vHit1mm.Y()+yt < (kSQ20_Y_size/2. + kThinYoffset)
				&& vHit1mm.Y()+yt > (-kSQ20_Y_size/2. + kThinYoffset)
		) {


			//coordinate in the thin detector plane
			xThin = vHit1mm.X()+xt;
			yThin = vHit1mm.Y()+yt;

			mapYbin = yThin*kSQL_Y_strips/kSQL_Y_size + kSQL_Y_strips/2.;	//ok
			mapXbin = xThin*kSQL_X_strips/kSQL_X_size - kSQL_X_strips/2.;	//ok


			//calculation of corrected energy
			//TODO: check the mapping here at second
			currThickness = hThickness->GetBinContent(mapXbin+1, mapYbin+1);

			for (Int_t x20 = 0; x20 < kSQL_20_strips; x20++) {
				if (SQ20E[x20] > 0/*SQ20timeMult==1*/ && currThickness<30.) {
					SQ20EcorrHit[x20] = (kSQ20_norm_thickness/currThickness)*SQ20E[x20];
					//when multiplicity larger than 1, we obtain incorrect information
					SQ20EcorrHitMult++;
//					cout << mapXbin << "\t" << mapYbin << "\t" << x20
//						<< "\t" << currThickness
//						<< "\t" << SQ20E[x20]  << "\t" << SQ20EcorrHit[x20] << "\t" << endl;

				}
			}//for*/

		}//if inside the map

		//end of MWPC wire multiplicity == 1

		//for MWPC cluster multiplicity == 1

		if (SQLXtimeMult==1 && SQLYtimeMult==1 && xtc>-50. && ytc>-50.) {
			vHit1mm.SetXYZ(x1mm-xtc, y1mm-ytc, z1mm);
			vHit1mm *= zThin/z1mm;
			angleLeft = vHit1mm.Angle(vNorm);
		}


		if (SQLXtimeMult==1 && SQLYtimeMult==1
				//check if the particle passed through sensitive area of thin detector
				&& vHit1mm.X()+xtc < (kSQ20_X_size/2. + kThinXoffset)
				&& vHit1mm.X()+xtc > (-kSQ20_X_size/2. + kThinXoffset)
				&& vHit1mm.Y()+ytc < (kSQ20_Y_size/2. + kThinYoffset)
				&& vHit1mm.Y()+ytc > (-kSQ20_Y_size/2. + kThinYoffset)
//				&& vHit1mm.X()+xtc>-30. && vHit1mm.X()+xtc<30.
//				&& vHit1mm.Y()+ytc>-30. && vHit1mm.Y()+ytc<30.
		) {


			xThin = vHit1mm.X()+xtc;
			yThin = vHit1mm.Y()+ytc;

			mapYbin = yThin*kSQL_Y_strips/kSQL_Y_size + kSQL_Y_strips/2.;	//ok
			mapXbin = xThin*kSQL_X_strips/kSQL_X_size - kSQL_X_strips/2.;	//ok
//			mapYbin = (vHit1mm.Y()+ytc+kThinYoffset+30.-0.5*60./16.)*16./60.;
//			mapXbin = (vHit1mm.X()+xtc+kThinXoffset+30.-0.5*60./32.)*32./60.;
			//coordinate in the thin detector plane


			//calculation of corrected energy
			currThickness = hThickness->GetBinContent(mapXbin+1, mapYbin+1);

			for (Int_t x20 = 0; x20 < kSQL_20_strips; x20++) {
				if (SQ20E[x20] > 0/*SQ20timeMult==1*/ && currThickness<30.) {
					SQ20EcorrHitC[x20] = (kSQ20_norm_thickness/currThickness)*SQ20E[x20];
					SQ20EcorrHitCMult++;
					//					cout << mapXbin << "\t" << mapYbin << "\t" << x20
					//						<< "\t" << currThickness
					//						<< "\t" << SQ20E[x20]  << "\t" << SQ20EcorrHit[x20] << "\t" << endl;

				}
			}//for*/

			if (SQ20EcorrHitCMult==1) {
				//TODO: check this statement
				expectedThinStrip = (xThin + 25. - kThinXoffset)*(16./50.)-0.5;
			}

		}//if inside the map

		//end of MWPC wire multiplicity == 1

		if (
				SQLXtimeMult>0 && SQLXEtimeFilteredSum>0
				&& trigger==3
				&& x1MultC==1 && y1MultC==1 && x2MultC==1 && y2MultC==1
		) {

			tw->Fill();
		}//if TTree::Fill

		if (eventNo%100000 == 0 && eventNo !=0) {
			cout << eventNo << " events processed." << endl;
		}
	}//for events

	cout << nevents << " entries processed." << endl;

	fw->cd();
	tw->Write();
	fw->Close();

}


void fillChain(const TString beam = "he", const Int_t from = 0, const Int_t to = 9, const Int_t noevents = 0) {

	TStopwatch stopwatch;
	stopwatch.Start();

//	cout << "aksjdlkajbdlkjasd" << endl;

	for (Int_t i = from; i <= to; i++) {
		fillTree(beam, i, noevents);
	}

	cout << "Finished in "<< stopwatch.RealTime() << " seconds" << endl;
	cout << "Finished in "<< stopwatch.RealTime()/60. << " minutes" << endl;
}