BeAnalysis.cpp

/*
 * BeAnalysis.cpp
 *
 *  Created on: Jul 19, 2017
 *      Author: vratik
 */

#include "BeAnalysis.h"

#include "BeWork.h"
#include "TCanvas.h"

BeAnalysis::BeAnalysis() : che(0), spectra(1) {
	// TODO Auto-generated constructor stub

	for (Int_t i = 0; i < 6; i++) {
		chs[i] = NULL;
		ti[i] = NULL;
	}

	SetCuts();
	SetChainsToDraw();

}

BeAnalysis::~BeAnalysis() {
	// TODO Auto-generated destructor stub
}

void BeAnalysis::OpenSimChains() {

	const Int_t lowSimFile[6] = {0, 0, 0, 0, 0, 0};
	const Int_t upSimFile[6] = {14, 14, 14, 14, 14, 14};

	TString chainAl0name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_al_0_50-85_";
	TString inputTreeAl0name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_al_0_50-85_000.root";
	//isotropic, 0 degrees
	TString chainNoAl0name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_noal_0_50-85_";
	TString inputTreeNoAl0name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_noal_0_50-85_000.root";
	//aligned, 180 degrees
	TString chainAl180name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_al_180_50-85_";
	TString inputTreeAl180name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_al_180_50-85_000.root";
	//isotropic, 180 degrees
	TString chainNoAl180name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_noal_180_50-85_";
	TString inputTreeNoAl180name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_noal_180_50-85_000.root";
	//isotropic, 90 degrees
	TString chainAl90name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_al_90_50-85_";
	TString inputTreeAl90name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_al_90_50-85_000.root";
	//isotropic, 90 degrees
	TString chainNoAl90name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_noal_90_50-85_";
	TString inputTreeNoAl90name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_noal_90_50-85_000.root";

//	TChain *chsAl0 = BeWork::OpenChain(chainAl0name.Data(),
	chs[0] = BeWork::OpenChain(chainAl0name.Data(),
			lowSimFile[0], upSimFile[0], "simbe", 4, "sbeam");
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing simulated data", chs[0]->GetEntries(), chs[0]->GetName());
	ti[0] = BeWork::OpenTree(inputTreeAl0name.Data(), "sbeam", 2);
	//	tiAl0->SetLineWidth(2);

//	TChain *chsNoAl0 = BeWork::OpenChain(chainNoAl0name.Data(),
	chs[1] = BeWork::OpenChain(chainNoAl0name.Data(),
			lowSimFile[1], upSimFile[1], "simbe", 4, "sbeam");
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing simulated data", chs[1]->GetEntries(), chs[1]->GetName());
	ti[1] = BeWork::OpenTree(inputTreeNoAl0name.Data(), "sbeam", 2);
	//	tiNoAl0->SetLineWidth(2);
	//
//	TChain *chsAl180 = BeWork::OpenChain(chainAl180name.Data(),
	chs[2] = BeWork::OpenChain(chainAl180name.Data(),
			lowSimFile[2], upSimFile[2], "simbe", 4, "sbeam");
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing simulated data", chs[2]->GetEntries(), chs[2]->GetName());
	ti[2] = BeWork::OpenTree(inputTreeAl180name.Data(), "sbeam", 2);
	//	tiAl180->SetLineWidth(2);

//	TChain *chsNoAl180 = BeWork::OpenChain(chainNoAl180name.Data(),
	chs[3] = BeWork::OpenChain(chainNoAl180name.Data(),
			lowSimFile[3], upSimFile[3], "simbe", 4, "sbeam");
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing simulated data", chs[3]->GetEntries(), chs[3]->GetName());
	ti[3] = BeWork::OpenTree(inputTreeNoAl180name.Data(), "sbeam", 2);
	//	tiNoAl180->SetLineWidth(2);

//	TChain *chsAl90 = BeWork::OpenChain(chainAl90name.Data(),
	chs[4] = BeWork::OpenChain(chainAl90name.Data(),
			lowSimFile[4], upSimFile[4], "simbe", 4, "sbeam");
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing simulated data", chs[4]->GetEntries(), chs[4]->GetName());
	ti[4] = BeWork::OpenTree(inputTreeAl90name.Data(), "sbeam", 2);
	//	tiAl90->SetLineWidth(2);

//	TChain *chsNoAl90 = BeWork::OpenChain(chainNoAl90name.Data(),
	chs[5] = BeWork::OpenChain(chainNoAl90name.Data(),
			lowSimFile[5], upSimFile[5], "simbe", 4, "sbeam");
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing simulated data", chs[5]->GetEntries(), chs[5]->GetName());
	ti[5] = BeWork::OpenTree(inputTreeNoAl90name.Data(), "sbeam", 2);
	//	tiNoAl90->SetLineWidth(2);

}

void BeAnalysis::OpenExpChain() {

	Int_t lowExpFile = 0;
	Int_t upExpFile = 40;

	//experimental chain
	che = BeWork::OpenChain("../../../be/rootdata/correlations/Be.",
			lowExpFile, upExpFile, "beonly");										//original file
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing experimental data",
			che->GetEntries(), che->GetName());

}

void BeAnalysis::SetCuts() {
	cBe20 = "fBeIM>0 && fBeIM<20";
	cBe3 = "fBeIM>0 && fBeIM<3";
	cBeWork = "fBeIM>0 && fBeIM<10";

	cBe0_14 = "fBeIM>0 && fBeIM<1.4";
	//	TCut cBe0_14 = "fBeIM>0 && fBeIM<1.6";
	//	TCut cBe0_14 = "fBeIM>0 && fBeIM<1.20";
	cBe14_19 = "fBeIM>1.4 && fBeIM<1.9";
	cBe19_25 = "fBeIM>1.9 && fBeIM<2.5";
	cBe25_31 = "fBeIM>2.5 && fBeIM<3.1";
	cBe31_37 = "fBeIM>3.1 && fBeIM<3.7";
	//array:
	cBeE[0] = cBe0_14;
	cBeE[1] = cBe14_19;
	cBeE[2] = cBe19_25;
	cBeE[3] = cBe25_31;
	cBeE[4] = cBe31_37;

	cEpsilonT = "fTpp/fBeIM<0.2";
	//	TCut cEpsilonT = "fTpp/fBeIM<1.";
	//	TCut cEpsilonY = "fTap/fBeIM<0.5";
	cEpsilonY = "fTap/fBeIM>0.7";

	//raw files
	//energy cuts
	crBe20 = "f6BeIM>0 && f6BeIM<20";
	crBe3 = "f6BeIM>0 && f6BeIM<3";
	crBeWork = "f6BeIM>0 && f6BeIM<10";

	crBe0_14 = "f6BeIM>0 && f6BeIM<1.4";
	//	TCut crBe0_14 = "f6BeIM>0 && f6BeIM<1.6";
	//	TCut crBe0_14 = "f6BeIM>0.5 && f6BeIM<1.20";
	crBe14_19 = "f6BeIM>1.4 && f6BeIM<1.9";
	crBe19_25 = "f6BeIM>1.9 && f6BeIM<2.5";
	crBe25_31 = "f6BeIM>2.5 && f6BeIM<3.1";
	crBe31_37 = "f6BeIM>3.1 && f6BeIM<3.7";
	//array:
	crBeE[0] = crBe0_14;
	crBeE[1] = crBe14_19;
	crBeE[2] = crBe19_25;
	crBeE[3] = crBe25_31;
	crBeE[4] = crBe31_37;
	//angular cuts


	//	TCut crEpsilonT = "fTpp/f6BeIM>0.8";
	crEpsilonT = "fTpp/f6BeIM<0.2";
	//	TCut crEpsilonT = "fTpp/f6BeIM<1.";
	//	TCut crEpsilonY = "fTap/f6BeIM<0.5";
	crEpsilonY = "fTap/f6BeIM>0.7";

	//simulation input
	ciBe0_14 = "E_IM>0 && E_IM<1.4";
	//	TCut ciBe0_14 = "E_IM>0 && E_IM<1.6";
	//	TCut ciBe0_14 = "E_IM>0.5 && E_IM<1.20";
	ciBe14_19 = "E_IM>1.4 && E_IM<1.9";
	ciBe19_25 = "E_IM>1.9 && E_IM<2.5";
	ciBe25_31 = "E_IM>2.5 && E_IM<3.1";
	ciBe31_37 = "E_IM>3.1 && E_IM<3.7";
	//array:
	ciBeE[5] = ciBe0_14;
	ciBeE[5] = ciBe14_19;
	ciBeE[5] = ciBe19_25;
	ciBeE[5] = ciBe25_31;
	ciBeE[5] = ciBe31_37;

	ciEpsilon = "sTpp/E_IM<0.2";
	//	TCut ciEpsilon = "sTpp/E_IM<1.";
	//	TCut ciEpsilonY = "sTap/E_IM<0.5";
	ciEpsilonY = "sTap/E_IM>0.7";

	cQ = "TMath::Abs(fQLiP)<10";
	cProtons = "fP1Lab.fE-938.272<50 && fP2Lab.fE-938.272<50";
	//	TCut cProtons = "fP1Lab.fE-938.272<34 && fP2Lab.fE-938.272<34 && fP1Lab.fE-938.272>1 && fP2Lab.fE-938.272>1";
}

void BeAnalysis::SetCMAngularRange(Int_t minAngle, Int_t maxAngle) {
	kMinAngle = minAngle;
	kMaxAngle = maxAngle;
}

void BeAnalysis::SetChainsToDraw(Bool_t ch0, Bool_t ch1, Bool_t ch2, Bool_t ch3, Bool_t ch4,
		Bool_t ch5) {

	kChains[0] = ch0;
	kChains[1] = ch1;
	kChains[2] = ch2;
	kChains[3] = ch3;
	kChains[4] = ch4;
	kChains[5] = ch5;

}

void BeAnalysis::SetRangeProportion(Float_t rangeProportion, Bool_t autoRange) {
	kRangeProportion = rangeProportion;
	kAutoRange = autoRange;
}

void BeAnalysis::SetNoExpEvents() {

	eMaxEvents = 5000000;

	eEventsAl0 = eMaxEvents;
	eEventsNoAl0 = eMaxEvents;
	eEventsAl180 = eMaxEvents;
	eEventsNoAl180 = eMaxEvents;
	eEventsAl90 = eMaxEvents;
	eEventsNoAl90 = eMaxEvents;

	eEvents[0] = eEventsAl0;
	eEvents[1] = eEventsNoAl0;
	eEvents[2] = eEventsAl180;
	eEvents[3] = eEventsNoAl180;
	eEvents[4] = eEventsAl90;
	eEvents[5] = eEventsNoAl90;

}

void BeAnalysis::SetNoSimEvents() {

	sMaxEvents = 3000000;

	sEventsAl0 = sMaxEvents;
	sEventsNoAl0 = 2700000;
	sEventsAl180 = sMaxEvents;
	sEventsNoAl180 = 2700000;
	sEventsAl90 = 2690000;
	sEventsNoAl90 = 2710000;

	sEvents[0] = sEventsAl0;
	sEvents[1] = sEventsNoAl0;
	sEvents[2] = sEventsAl180;
	sEvents[3] = sEventsNoAl180;
	sEvents[4] = sEventsAl90;
	sEvents[5] = sEventsNoAl90;

}

void BeAnalysis::SetRatiosGStoEX() {

	sRatioAl0 = "sRatio>0.051 && sRatio<0.101";
	//	TCut sRatioAl0 = "sRatio>0.060 && sRatio<0.110";
	sRatioNoAl0 = "sRatio>0.050 && sRatio<0.100";
	sRatioAl180 = "sRatio>0.047 && sRatio<0.097";
	sRatioNoAl180 = "sRatio>0.048 && sRatio<0.098";
	sRatioAl90 = "sRatio>0.051 && sRatio<0.101";
	sRatioNoAl90 = "sRatio>0.051 && sRatio<0.101";

	sRatio[0] = sRatioAl0;
	sRatio[1] = sRatioNoAl0;
	sRatio[2] = sRatioAl180;
	sRatio[3] = sRatioNoAl180;
	sRatio[4] = sRatioAl90;
	sRatio[5] = sRatioNoAl90;

}

void BeAnalysis::SetCMAngularCuts() {
	cAngles = "fBeThetaCM1>60*TMath::DegToRad() && fBeThetaCM1<75*TMath::DegToRad()";
	crAngles = "f6BeThetaCM1>60*TMath::DegToRad() && f6BeThetaCM1<75*TMath::DegToRad()";
}

void BeAnalysis::Spectra() {

	if (!spectra) {
		return;
	}

	if (!che) return;

	TCanvas *cSpectra = new TCanvas();
	TString canvasTitle;

	TH1F *hsSpectra[6];
	TH1F *heSpectra[6];
	TH1F *hSdiff[6];

	canvasTitle.Form("energy spectra fitting;\t(%d,%d) degrees", kMinAngle, kMaxAngle);
	cSpectra->SetTitle(canvasTitle.Data());
	cSpectra->Divide(2, 3);

	che->SetLineColor(kGreen);

	for (Int_t j = 0; j < 6; j++) {	//different files
		if (!kChains[j]) continue;
		cout << kChains[j] << endl;
		cSpectra->cd(j+1);
		chs[j]->SetLineColor(kBlue);
		chs[j]->SetFillColor(kWhite);
		hsName.Form("hsSpectra%d", j);
		drawCommand.Form("f6BeIM>>%s(200,0,10)", hsName.Data());
		//			chs[j]->Draw(drawCommand.Data(), cQ && crBeWork && crAngles && sRatio[j], "", sEvents[j]);
		chs[j]->Draw(drawCommand.Data(), cProtons && cQ && crBeWork && crAngles && sRatio[j], "", sEvents[j]);
		hsSpectra[j] = (TH1F*)gPad->FindObject(hsName.Data());
		heName.Form("heSpectra%d", j);
		drawCommand.Form("fBeIM>>%s(200,0,10)", heName.Data());
		//			che->Draw(drawCommand.Data(), cQ && cBeWork && cAngles, "same", eEvents[j]);
		che->Draw(drawCommand.Data(), cProtons && cQ && cBeWork && cAngles, "", eEvents[j]);
		heSpectra[j] = (TH1F*)gPad->FindObject(heName.Data());

		//			heSpectra[j]->Draw("");
		heSpectra[j]->Draw("E");
		heSpectra[j]->SetLineColor(1);
		heSpectra[j]->SetTitle("");
		heSpectra[j]->SetXTitle("E_{T} [MeV]");
		heSpectra[j]->GetXaxis()->SetTitleOffset(0.95);
		heSpectra[j]->GetXaxis()->CenterTitle();
		heSpectra[j]->SetYTitle("counts");
		heSpectra[j]->GetYaxis()->SetTitleOffset(0.7);
		heSpectra[j]->GetYaxis()->CenterTitle();
		hsSpectra[j]->SetLineColor(kGray+1);
		hsSpectra[j]->SetFillColor(kGray+1);
		hsSpectra[j]->Draw("same");
		heSpectra[j]->Draw("same E");

		if (kAutoRange) {
			Float_t leftMaxMC = kRangeProportion*hsSpectra[j]->GetMaximum();
			Float_t leftMaxE = kRangeProportion*heSpectra[j]->GetMaximum();
			//				hsSpectra[j]->GetYaxis()->SetRangeUser(0, leftMaxE);
			leftMaxMC > leftMaxE ? heSpectra[j]->GetYaxis()->SetRangeUser(0, leftMaxMC) : heSpectra[j]->GetYaxis()->SetRangeUser(0, leftMaxE);
		}


		hSdiff[j] = new TH1F(*heSpectra[j]);
		hSdiff[j]->Add(hsSpectra[j], -1);
		hSdiff[j]->SetLineColor(kRed);
		//			hSdiff[j]->Draw("same");



		cSpectra->Update();

		cout << endl;
		Info("statesRatioFitting", "MC %d", j+1);
		//			Info("statesRatioFitting", "ground state: %3.1f (exp)/ %3.1f (sim) = %3.3f", heSpectra[j]->Integral(0,20), hsSpectra[j]->Integral(0,20), heSpectra[j]->Integral(0,20)/hsSpectra[j]->Integral(0,20) );
		//			Info("statesRatioFitting", "left slope: %3.1f (exp)/ %3.1f (sim) = %3.3f\n", heSpectra[j]->Integral(20,32), hsSpectra[j]->Integral(20,32), heSpectra[j]->Integral(20,32)/hsSpectra[j]->Integral(20,32) );
		//0-2 MeV, 2-3.1 MeV
		//			Info("statesRatioFitting", "ground state: %3.1f (exp)/ %3.1f (sim) = %3.3f", heSpectra[j]->Integral(0,40), hsSpectra[j]->Integral(0,40), heSpectra[j]->Integral(0,40)/hsSpectra[j]->Integral(0,40) );
		//			Info("statesRatioFitting", "left slope: %3.1f (exp)/ %3.1f (sim) = %3.3f\n", heSpectra[j]->Integral(40,64), hsSpectra[j]->Integral(40,64), heSpectra[j]->Integral(40,64)/hsSpectra[j]->Integral(40,64) );
		//0-2 MeV, 2.5-3.5 MeV
		Info("statesRatioFitting", "ground state: %3.1f (sim)/ %3.1f (exp) = %3.3f", hsSpectra[j]->Integral(0,40), heSpectra[j]->Integral(0,40), hsSpectra[j]->Integral(0,40)/heSpectra[j]->Integral(0,40) );
		Info("statesRatioFitting", "left slope: %3.1f (sim)/ %3.1f (exp) = %3.3f\n", hsSpectra[j]->Integral(50,70), heSpectra[j]->Integral(50,70), hsSpectra[j]->Integral(50,70)/heSpectra[j]->Integral(50,70) );
		//0-1.3 MeV, 2-3.1 MeV
		//			Info("statesRatioFitting", "ground state: %3.1f (exp)/ %3.1f (sim) = %3.3f", heSpectra[j]->Integral(0,26), hsSpectra[j]->Integral(0,26), heSpectra[j]->Integral(0,26)/hsSpectra[j]->Integral(0,26) );
		//			Info("statesRatioFitting", "left slope: %3.1f (exp)/ %3.1f (sim) = %3.3f\n", heSpectra[j]->Integral(40,64), hsSpectra[j]->Integral(40,64), heSpectra[j]->Integral(40,64)/hsSpectra[j]->Integral(40,64) );
	}//for j

	/*if (savePictures) {
			canvasName.Form("%sfigSpectra:%s%s", ppath.Data(), configuration.Data(), pictFormat.Data());
			cout << canvasName.Data() << endl;
			cSpectra->SaveAs(canvasName.Data());
			cSpectra->cd();
			cSpectra->Close();
		}//if save pictures*/


	Info("sfAngInt_spectra.cxx", "Finished.");

}