diff --git a/Be/BeAnalysis.cpp b/Be/BeAnalysis.cpp index 7866354a512ee055f9fd94cf3d15f59e703c7288..9fc7353fe899da35e94c18a02ca5c42e70a320ce 100644 --- a/Be/BeAnalysis.cpp +++ b/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."); + +} diff --git a/Be/BeAnalysis.h b/Be/BeAnalysis.h index ee610520d90a7a9ab95c2aced9bf57acf952147c..4421ef8815d74a8995a846858a5edc6a533c3943 100644 --- a/Be/BeAnalysis.h +++ b/Be/BeAnalysis.h @@ -12,6 +12,7 @@ //#include #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_ */ diff --git a/macros/BeCorrPRC/lib_test.cxx b/macros/BeCorrPRC/lib_test.cxx new file mode 100755 index 0000000000000000000000000000000000000000..0d7b2a5bf6dc98d452b778a88677b0b27d8a07f9 --- /dev/null +++ b/macros/BeCorrPRC/lib_test.cxx @@ -0,0 +1,15 @@ +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(); +}