fillChain.cxx

#include "TSystem.h"
#include "TFile.h"
#include "TTree.h"
#include "TVector3.h"
//#include "../src/TNeEvent.h"
//#include "../../AculUtils/TELoss/TELoss.h"

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

void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {

	TString inFile;
//	inFile.Form("~/data/exp1804/h5_14_%02d.root", nofile);
	inFile.Form("~/data/exp1804/be10_03_%d0.root", nofile);			//files 00,10,...,60
//	inFile.Form("~/data/exp1804/be10_%02d_00.root", nofile);		//files 01,02, 03, 05
//	inFile.Form("~/data/exp1804/calib/si_20_03.root");

	TString outFile;
//	outFile.Form("~/data/exp1804/h5_14_%02d_calib.root", nofile);
	outFile.Form("~/data/exp1804/be10_03_%d0_calib.root", nofile);	//files 00,10,...,60
//	outFile.Form("~/data/exp1804/be10_%02d_00_calib.root", nofile); //files 01,02, 03, 05
//	outFile.Form("~/data/exp1804/calib/si_20_03_calib.root");

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

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

	TNeEvent *revent = new TNeEvent();

	tr->SetBranchAddress("NeEvent.", &revent);

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

	Float_t SQ20E[16];
	Float_t SQ20Ecorr[16];
	Float_t SQ20EcorrHit[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;

	Float_t SQRXE[32];
	Float_t SQRXEsum;
	Int_t SQRXmult;
	Float_t TOF, dEbeam;
	Float_t x1p, y1p, x2p, y2p, xt, yt;


	//left 1 mm position
	Float_t x1mm, y1mm;
	const Float_t z1mm = 230.;

	Float_t xThin, yThin;
	const Float_t zThin = 230.-53.6;

	const Float_t xThinOffset = -3., yThinOffset = -1.8;

	Int_t mapXbin, mapYbin;

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

	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("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("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("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");


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

//	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();

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

	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");

//	std::cout << std::setprecision(1) << std::fixed;


	const Int_t kSQL_X_strips = 32;
	const Int_t kSQL_Y_strips = 16;
	const Int_t kSQL_20_strips = 16;

	const Double_t kSQLY_energy_thr = 1.;
	const Double_t kSQLX_energy_thr = 1.;
	const Double_t kSQL20_energy_thr = 1.2;

	const Double_t kSQ20_norm_thickness = 20.;

//	for (Int_t yi = 0; yi < kSQL_Y_strips; yi++) {
//		cout << "y bin: " << yi+1 << "\t\t";
//		for (Int_t xi = 0; xi < kSQL_X_strips; xi++) {
//			cout << hThickness->GetBinContent(xi+1, yi+1) << "\t";
//			if (xi == kSQL_X_strips-1) cout << endl;
//		}
//	}

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

	fw->cd();

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

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

//		cout << eventNo << endl;

		tr->GetEvent(eventNo);

//		cout << eventNo << endl;
//		cout << revent->SQX_L[0] << endl;
//		SQLXE[0] = revent->SQX_L[0];

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

		SQLXEtimeFilteredSum = 0.;
		SQLYEtimeFilteredSum = 0.;

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

		SQ20timeMult = 0;
		SQLXtimeMult = 0;
//		SQLXtimeSum = 0;
		SQLYtimeMult = 0;



//cout << endl;

		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.;

			energy = pSQX_L_EC.Energy(revent->SQX_L[i], i);

//			if (i == 5 && i == 6) continue;

			if (energy > kSQLX_energy_thr) {
				SQLXE[i] = energy;
				SQLXEsum += SQLXE[i];
				SQLXmult++;
			}

			if (i<16 && i!=5 && i!=6 && i!=0 && i!=1
					&& SQLXE[i]>kSQLX_energy_thr
					&& revent->tSQX_L[i]*0.3>320 && revent->tSQX_L[i]*0.3<380) {
				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]>kSQLX_energy_thr
					&& revent->tSQX_L[i]*0.3>410 && revent->tSQX_L[i]*0.3<470) {
				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 (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.;
			//				SQRXE[i] = 0;

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

			energy = pSQ20_EC.Energy(revent->SQ20[i], i);

//			if (energy > 1.) {
				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]>kSQL20_energy_thr && SQ20time[i]>400 && SQ20time[i]<500) SQ20timeMult++;


			if (i==0) continue;

			energy = pSQY_L_EC.Energy(revent->SQY_L[i], i);

			//				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 (energy > kSQLY_energy_thr) {
				SQLYE[i] = energy;
				SQLYEsum += SQLYE[i];
				SQLYmult++;
			}

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

			if (SQLYE[i]>kSQLY_energy_thr && time>325 && time<333) {
//				cout << "!!!!!!!!!!!!" << endl;
				SQLYtime[i] = time;
				SQLYtimeMult++;
				SQLYEtimeFiltered[i]=energy;
				SQLYEtimeFilteredSum = SQLYEtimeFilteredSum + SQLYEtimeFiltered[i];
			}



			//				energy = pSQX_R_EC.Energy(revent->SQX_R[i], i);
			//
			//				if (energy > 1.0) {
			//					SQRXE[i] = energy;
			//					SQRXEsum += SQRXE[i];
			//					SQRXmult++;
			//				}
		}//for 16

//		cout << SQLXtimeMult << endl;

//		cout << endl;

		///////////////////////////////////////////////
		//correction for thickness inhomogeneity
		///////////////////////////////////////////////

		Double_t currThickness;

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

			for (Int_t yi = 0; yi < kSQL_Y_strips; yi++) {
				//			cout << "y bin: " << yi+1 << "\t\t";
				if (SQLYE[yi]>kSQLY_energy_thr) {

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

						if (SQLXE[xi] > kSQLX_energy_thr) {
							currThickness = hThickness->GetBinContent(xi+1, yi+1);
//							if (currThickness<30.) {
//								Int_t probableThinStrip = xi-2/2;
//								SQ20Ecorr[0] =
//							}

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

//									currThickness = kSQ20_norm_thickness + (currThickness - kSQ20_norm_thickness)*0.5;

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

//									cout << xi << "\t" << yi << "\t" << x20
//											<< "\t" << currThickness
//											<< "\t" << SQ20E[x20]  << "\t" << SQ20Ecorr[x20] << "\t" << endl;
								}
							}//for*/
						}



						//					cout << hThickness->GetBinContent(xi+1, yi+1) << "\t";
						//					if (xi == kSQL_X_strips-1) cout << endl;
					}
				}
			}

		}//if for correction according energy multiplicity


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

		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;
		}

		const Float_t l12 = 546.;
		const Float_t lt = 270.;

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

		if (revent->x1[0]<1000 && revent->y1[0]<1000 && revent->nx1==1 && revent->ny1==1) {
			x1p = revent->x1[0]*1.25-20.;
			y1p = revent->y1[0]*1.25-20.;
		}

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

		if (x1p > -50. && y1p > -50. && x2p > -50. && y2p > -50.) {
			xt = x1p - (x1p - x2p)*(l12/lt);
			yt = y1p - (y1p - y2p)*(l12/lt);
		}



		///////////////////////////////////////////////
		//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.;
			}
		}
//		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.;
			}
		}


		if (SQLXtimeMult==1 && SQLYtimeMult==1 && xt>-50. && yt>-50.) {
			vHit1mm.SetXYZ(x1mm-xt, y1mm-yt, z1mm);
//			cout << vHit1mm.X() << "\t" << vHit1mm.Y() << endl;
//			cout << "\t" << zThin/z1mm*vHit1mm.X() << "\t" << zThin/z1mm*vHit1mm.Y() << endl;
			vHit1mm *= zThin/z1mm;
//			cout << "\t" << vHit1mm.X() << "\t" << vHit1mm.Y() << "\t" << vHit1mm.Z() << endl;
//			cout << "\t" << xt << "\t" << yt << "\t" << 0 << endl;
//			cout << zThin/z1mm << endl;
			angleLeft = vHit1mm.Angle(vNorm);
//			xThin = vHit1mm.X()+xt;
//			yThin = vHit1mm.Y()+yt;

//			cout << angleLeft << endl;
		}

//		cout << angleLeft << endl;

//		cout << x1mm << "\t" << y1mm << endl;

//		cout << vHit1mm.X() << "\t" << vHit1mm.Y() << endl;



		if (SQLXtimeMult==1 && SQLYtimeMult==1
				&& vHit1mm.X()+xt>-30. && vHit1mm.X()+xt<30.
				&& vHit1mm.Y()+yt>-30. && vHit1mm.Y()+yt<30.) {
			mapYbin = (vHit1mm.Y()+yt+30.-0.5*60./16.)*16./60.;
			mapXbin = (vHit1mm.X()+xt+30.-0.5*60./32.)*32./60.;
			xThin = vHit1mm.X()+xt;
			yThin = vHit1mm.Y()+yt;
//			cout << "\t\t" << xThin << "\t" << yThin << endl;

//			cout << mapYbin << "\t" << "\t" << y1mm << "\t" << xt << endl;

			//calculation of corrected energy
			currThickness = hThickness->GetBinContent(mapXbin+1, mapYbin+1);
			for (Int_t x20 = 0; x20 < kSQL_20_strips; x20++) {
				if (SQ20E[x20] > kSQL20_energy_thr/*SQ20timeMult==1*/ && currThickness<30.) {
					SQ20EcorrHit[x20] = (kSQ20_norm_thickness/currThickness)*SQ20E[x20];
					SQ20EcorrHitMult++;
//					cout << mapXbin << "\t" << mapYbin << "\t" << x20
//						<< "\t" << currThickness
//						<< "\t" << SQ20E[x20]  << "\t" << SQ20EcorrHit[x20] << "\t" << endl;

				}
			}//for*/

		}//if inside the map

//		cout << mapYbin << "\t" << "\t" << y1mm << endl;
//		cout << mapXbin << "\t" << xThin << "\t" << x1mm << endl;


//		cout << mapXbin << "\t" << mapYbin << endl;

		tw->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 Int_t from = 0, const Int_t to = 9, const Int_t noevents = 0) {
	for (Int_t i = from; i <= to; i++) {
		fillTree(i, noevents);
	}
}