BeAnalysis::Spectra() implemented

Be/BeAnalysis.cpp

@@ -6,11 +6,21 @@
  */
 
 #include "BeAnalysis.h"
+
 #include "BeWork.h"
+#include "TCanvas.h"
 
-BeAnalysis::BeAnalysis() {
+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() {
@@ -89,11 +99,271 @@ void BeAnalysis::OpenExpChain() {
 	Int_t lowExpFile = 0;
 	Int_t upExpFile = 40;
 
-//	cout << "aklsdjkajhsdkajshd" << endl;
 	//experimental chain
-	TChain *che = BeWork::OpenChain("../../../be/rootdata/correlations/Be.",
+	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.");
+
+}

Be/BeAnalysis.h

@@ -12,6 +12,7 @@
 
 //#include <TObject.h>
 #include "TChain.h"
+#include "TCut.h"
 
 using std::cout;
 using std::endl;
@@ -28,11 +29,128 @@ public:
 
 	TChain* GetExpChain() {return che;};
 	TChain* GetSimChain(Int_t i) {return chs[i];};
+	void SetChainsToDraw(Bool_t ch0 = 1, Bool_t ch1 = 1, Bool_t ch2 = 1, Bool_t ch3 = 1, Bool_t ch4 = 1, Bool_t ch5 = 1);
+
+	void SetCuts();
+
+	void SetCMAngularRange(Int_t minAngle, Int_t maxAngle);
+	void SetCMAngularCuts();
+	void SetRangeProportion(Float_t rangeProportion = 1.1, Bool_t autoRange = 1);
+
+	void SetNoExpEvents();
+	void SetNoSimEvents();
+	void SetRatiosGStoEX();
+
+	void Spectra();
 //
 private:
 	TChain *che;		//chain with experimental trees
 	TChain *chs[6];		//array of chains with simulation trees
 	TTree *ti[6];		//array of chains with simulation input
+
+	Int_t kMinAngle;
+	Int_t kMaxAngle;
+
+	TCut cAngles;
+	TCut crAngles;
+
+	Bool_t kChains[6];
+
+	Bool_t kAutoRange;
+	Float_t kRangeProportion;
+
+	Long64_t eMaxEvents;
+	Long64_t eEventsAl0;
+	Long64_t eEventsNoAl0;
+	Long64_t eEventsAl180;
+	Long64_t eEventsNoAl180;
+	Long64_t eEventsAl90;
+	Long64_t eEventsNoAl90;
+	Long64_t eEvents[6];
+
+	Long64_t sMaxEvents;
+	Long64_t sEventsAl0;
+	Long64_t sEventsNoAl0;
+	Long64_t sEventsAl180;
+	Long64_t sEventsNoAl180;
+	Long64_t sEventsAl90;
+	Long64_t sEventsNoAl90;
+	Long64_t sEvents[6];
+
+	TCut sRatioAl0;
+	TCut sRatioNoAl0;
+	TCut sRatioAl180;
+	TCut sRatioNoAl180;
+	TCut sRatioAl90;
+	TCut sRatioNoAl90;
+	TCut sRatio[6];
+
+	//cuts
+	TCut cBe20;
+	TCut cBe3;
+	TCut cBeWork;
+	TCut cBe0_14;
+	//	TCut cBe0_14 = "fBeIM>0 && fBeIM<1.6";
+	//	TCut cBe0_14 = "fBeIM>0 && fBeIM<1.20";
+	TCut cBe14_19;
+	TCut cBe19_25;
+	TCut cBe25_31;
+	TCut cBe31_37;
+	TCut cBeE[5];
+
+	TCut cEpsilonT;
+	//	TCut cEpsilonT = "fTpp/fBeIM<1.";
+	//	TCut cEpsilonY = "fTap/fBeIM<0.5";
+	TCut cEpsilonY;
+
+	//raw files
+	//energy cuts
+	TCut crBe20;
+	TCut crBe3;
+	TCut crBeWork;
+	TCut crBe0_14;
+	//	TCut crBe0_14 = "f6BeIM>0 && f6BeIM<1.6";
+	//	TCut crBe0_14 = "f6BeIM>0.5 && f6BeIM<1.20";
+	TCut crBe14_19;
+	TCut crBe19_25;
+	TCut crBe25_31;
+	TCut crBe31_37;
+	TCut crBeE[5];
+	//angular cuts
+
+
+	//	TCut crEpsilonT = "fTpp/f6BeIM>0.8";
+	TCut crEpsilonT;
+	//	TCut crEpsilonT = "fTpp/f6BeIM<1.";
+	//	TCut crEpsilonY = "fTap/f6BeIM<0.5";
+	TCut crEpsilonY;
+
+	//simulation input
+	TCut ciBe0_14;
+	//	TCut ciBe0_14 = "E_IM>0 && E_IM<1.6";
+	//	TCut ciBe0_14 = "E_IM>0.5 && E_IM<1.20";
+	TCut ciBe14_19;
+	TCut ciBe19_25;
+	TCut ciBe25_31;
+	TCut ciBe31_37;
+	TCut ciBeE[5];
+
+	TCut ciEpsilon;
+	//	TCut ciEpsilon = "sTpp/E_IM<1.";
+	//	TCut ciEpsilonY = "sTap/E_IM<0.5";
+	TCut ciEpsilonY;
+
+	TCut cQ;
+	TCut cProtons;
+	//	TCut cProtons = "fP1Lab.fE-938.272<34 && fP2Lab.fE-938.272<34 && fP1Lab.fE-938.272>1 && fP2Lab.fE-938.272>1";
+
+	Bool_t spectra;
+
+	//auxiliary strings
+	TString drawCommand;
+	TString hsName;
+	TString heName;
+
 };
 
 #endif /* BE_BEANALYSIS_H_ */

macros/BeCorrPRC/lib_test.cxx

+void lib_test() {
+	BeAnalysis ana;
+	ana.OpenExpChain();
+	ana.OpenSimChains();
+
+	ana.SetCMAngularRange(60, 75);
+	ana.SetCMAngularCuts();
+
+	ana.SetRangeProportion();
+	ana.SetNoExpEvents();
+	ana.SetNoSimEvents();
+	ana.SetRatiosGStoEX();
+
+	ana.Spectra();
+}