From 2abfe7e04ad748befa41640851a0b736cbdd1a3a Mon Sep 17 00:00:00 2001
From: morfouace <pierre.morfouace@cea.fr>
Date: Tue, 23 Jul 2024 17:12:28 +0200
Subject: [PATCH] updating analysis and macro for e850
---
Projects/Inelastic/Analysis.cxx | 2 +-
Projects/ana_e850/Analysis.cxx | 9 +-
Projects/ana_e850/macro/Mass/FillHistoMass.C | 10 +-
Projects/ana_e850/macro/Mass/FitMass.C | 81 ++++----
Projects/ana_e850/macro/Mass/GetTotalYield.C | 143 ++++++++++++++
.../ana_e850/macro/Mass/GetYieldEvolution.C | 179 ++++++++++++++++++
.../ana_e850/macro/Mass/PlotYieldEvolution.C | 20 +-
Projects/ana_e850/macro/Mass/check_gamma.C | 2 +
8 files changed, 402 insertions(+), 44 deletions(-)
create mode 100644 Projects/ana_e850/macro/Mass/GetTotalYield.C
create mode 100644 Projects/ana_e850/macro/Mass/GetYieldEvolution.C
diff --git a/Projects/Inelastic/Analysis.cxx b/Projects/Inelastic/Analysis.cxx
index 452e686e3..8ddfc1bbc 100644
--- a/Projects/Inelastic/Analysis.cxx
+++ b/Projects/Inelastic/Analysis.cxx
@@ -44,7 +44,7 @@ void Analysis::Init(){
InitInputBranch();
InitOutputBranch();
- my_Reaction = new NPL::Reaction("1n(238U,1n)238U@2.1");
+ my_Reaction = new NPL::Reaction("1n(56Fe,1n)56Fe@1.1");
neutron = new NPL::Nucleus("1n");
}
diff --git a/Projects/ana_e850/Analysis.cxx b/Projects/ana_e850/Analysis.cxx
index 568f26411..0aa090f97 100644
--- a/Projects/ana_e850/Analysis.cxx
+++ b/Projects/ana_e850/Analysis.cxx
@@ -398,14 +398,15 @@ void Analysis::VamosAnalysis(){
Tracking->CalculateReconstruction(FPMW->Xf, 1000*FPMW->Thetaf, m_Brho_ref, FF_Brho, Theta, FF_Path);
// FF_Path is in cm !
+ FF_Y1 = FPMW->PositionY[0];
+ FF_Y3 = FPMW->PositionY[2];
+
// T13 //
double path1 = FPMW->GetDetectorPositionZ(0)/10./cos(FPMW->Theta_in)/cos(FPMW->Phi_in);
double path2 = (FPMW->GetDetectorPositionZ(2)-7600)/10./cos(FPMW->Thetaf);
- FF_Y1 = FPMW->PositionY[0];
- FF_Y3 = FPMW->PositionY[2];
- FF_D13 = FF_Path - path1 + path2;
+ FF_D13 = FF_Path - path1 + path2;
FF_T13 = T13 - 20 + m_T13_Offset[FPMWPat];
- FF_V13 = FF_D13/FF_T13;
+ FF_V13 = (FF_D13/FF_T13)*(1-3e-5*FF_Y3);
FF_Beta13 = FF_V13/29.9792458;
FF_Gamma13 = 1./sqrt(1.0 - FF_Beta13*FF_Beta13);
FF_AoQ13 = 1.0*(FF_Brho/3.10761/FF_Beta13/FF_Gamma13);
diff --git a/Projects/ana_e850/macro/Mass/FillHistoMass.C b/Projects/ana_e850/macro/Mass/FillHistoMass.C
index ab8e9d206..2fb7782a7 100644
--- a/Projects/ana_e850/macro/Mass/FillHistoMass.C
+++ b/Projects/ana_e850/macro/Mass/FillHistoMass.C
@@ -62,6 +62,7 @@ void FillHistoMass(string part="part1"){
double VAMOS_TS_hour;
int Telescope;
int m_2alpha;
+ int FPMW_Section;
chain->SetBranchStatus("DeltaEcorr","true");
chain->SetBranchAddress("DeltaEcorr",&DeltaEcorr);
@@ -77,6 +78,9 @@ void FillHistoMass(string part="part1"){
chain->SetBranchAddress("Telescope",&Telescope);
chain->SetBranchStatus("m_2alpha","true");
chain->SetBranchAddress("m_2alpha",&m_2alpha);
+ chain->SetBranchStatus("FPMW_Section","true");
+ chain->SetBranchAddress("FPMW_Section",&FPMW_Section);
+
TH1F *hmass[11];
TH2F *hPID[8];
@@ -91,7 +95,7 @@ void FillHistoMass(string part="part1"){
//TH2F* hPID = new TH2F("PID","PID",1000,0,200,1000,0,60);
TH1F* hEx = new TH1F("Ex","Ex",500,0,20);
TH1F* hMassAll = new TH1F("hmass_all","hmass_all",1000,80,160);
-
+ TH2D* hMassVsEx = new TH2D("hmass_vs_ex","hmass_vs_ex",20,0,20,1000,80,160);
int nentries = chain->GetEntries();
cout << "**** Number of entries to be treated: " << nentries << " ****" << endl;
@@ -106,9 +110,10 @@ void FillHistoMass(string part="part1"){
if(Telescope>0) {
if(m_2alpha==0) hPID[Telescope-1]->Fill(Elab,DeltaEcorr);
- if(cut10Be[Telescope-1]->IsInside(Elab,DeltaEcorr)){
+ if(FPMW_Section>4 && cut10Be[Telescope-1]->IsInside(Elab,DeltaEcorr)){
hEx->Fill(Ex240Pu);
hMassAll->Fill(Mass);
+ hMassVsEx->Fill(Ex240Pu,Mass);
if(Ex240Pu>5.5 && Ex240Pu<6.5) hmass[0]->Fill(Mass);
else if(Ex240Pu>6.5 && Ex240Pu<7.5) hmass[1]->Fill(Mass);
else if(Ex240Pu>7.5 && Ex240Pu<8.5) hmass[2]->Fill(Mass);
@@ -132,6 +137,7 @@ void FillHistoMass(string part="part1"){
}
hEx->Write();
hMassAll->Write();
+ hMassVsEx->Write();
for(int i=0; i<8; i++)
hPID[i]->Write();
ofile->Close();
diff --git a/Projects/ana_e850/macro/Mass/FitMass.C b/Projects/ana_e850/macro/Mass/FitMass.C
index 62513072d..8df3169fa 100644
--- a/Projects/ana_e850/macro/Mass/FitMass.C
+++ b/Projects/ana_e850/macro/Mass/FitMass.C
@@ -1,18 +1,21 @@
void Fit(TH1F* hmass, string Energy);
TFile* ofile;
+int rebin=1;
//////////////////////////////////////////////////////
void FitMass(){
TFile* ifile = new TFile("histo_mass_all.root","read");
- //TFile* ifile = new TFile("histo_mass_orig.root","read");
//TFile* ifile = new TFile("histo_mass_part2.root","read");
ofile = new TFile("histo_mass_fitted.root","recreate");
- for(int i=0; i<10; i++){
+ for(int i=0; i<11; i++){
//for(int i=2; i<3; i++){
int Ex = i+6;
+ if(i>4) rebin=2;
+ else rebin=1;
+
string Energy = to_string(Ex) + "MeV";
TString histo_name = Form("hmass_%iMeV",i+6);
TH1F* h1 = (TH1F*)ifile->FindObjectAny(histo_name);
@@ -24,12 +27,14 @@ void FitMass(){
//////////////////////////////////////////////////////
void Fit(TH1F* hmass,string Energy){
- hmass->Rebin(2);
- hmass->Draw();
+ hmass->Rebin(rebin);
+ hmass->Draw("E1");
- TGraph* gyield = new TGraph();
int Amin = 82;
int Amax = 153;
+ //int Amin = 95;
+ //int Amax = 106;
+
int A = Amin;
@@ -45,27 +50,27 @@ void Fit(TH1F* hmass,string Energy){
for(int i=0; i<NumberOfA; i++){
g[i] = new TF1(Form("g%i",i),"gaus",A-0.4,A+0.4);
if(i==0){
- g[i]->SetParLimits(0,0.1,20000);
- g[i]->SetParLimits(1,A-0.2,A+0.2);
- g[i]->SetParLimits(2,0.1,0.4);
+ g[i]->SetParameter(0,100);
+ g[i]->SetParameter(1,Amin);
+ g[i]->SetParameter(2,0.3);
}
else if(i>0){
g[i]->SetParameter(0,g[i-1]->GetParameter(0));
- g[i]->SetParameter(1,g[i-1]->GetParameter(1));
+ g[i]->SetParameter(1,A);
g[i]->SetParameter(2,g[i-1]->GetParameter(2));
-
- g[i]->SetParLimits(1,A-0.2,A+0.2);
- g[i]->SetParLimits(2,0.1,0.3);
-
}
-
- hmass->Fit(g[i],"qR");
+
+ g[i]->SetParLimits(1,A-0.1,A+0.1);
+ g[i]->SetParLimits(2,0.25,0.4);
+
+
+ hmass->Fit(g[i],"qBR");
g[i]->Draw("lsame");
g[i]->GetParameters(¶[3*i]);
+ //cout << A << " " << g[i]->GetParameter(2) << endl;
A++;
-
}
@@ -76,38 +81,48 @@ void Fit(TH1F* hmass,string Energy){
}
TF1* total = new TF1("total",total_func,Amin-0.5, Amin+NumberOfA-0.5);
-
+ total->SetParameters(para);
A = Amin;
for(int i=0; i<NumberOfA; i++){
total->SetParameter(3*i,g[i]->GetParameter(0));
total->SetParameter(3*i+1,g[i]->GetParameter(1));
total->SetParameter(3*i+2,g[i]->GetParameter(2));
total->SetParLimits(3*i+1,A-0.2,A+0.2);
- total->SetParLimits(3*i+2,0.1,0.4);
- //cout << A << " " << g[i]->GetParameter(2) << endl;
+ total->SetParLimits(3*i+2,g[i]->GetParameter(2)*0.8,g[i]->GetParameter(2)*1.1);
A++;
}
- total->SetParameters(para);
total->SetLineColor(4);
- hmass->Fit(total,"Rq");
+ hmass->Fit(total,"RBq");
A = Amin;
+ //TF1* one_gauss = new TF1("one_gaus","gaus");
for(int i=0; i<NumberOfA; i++){
double Amplitude = total->GetParameter(3*i);
double mean = total->GetParameter(3*i+1);
double sigma = total->GetParameter(3*i+2);
double Amplitude_err = total->GetParError(3*i);
+ double mean_err = total->GetParError(3*i+1);
double sigma_err = total->GetParError(3*i+2);
- /*cout << "A= " << A << endl;
- cout << "Amplitude= " << Amplitude << endl;
- cout << "Mean= " << mean << endl;
- cout << "Sigma= " << sigma << endl;
-*/
- Integral[i] = Amplitude*sigma*sqrt(2*TMath::Pi());
- Integral_err[i] = sqrt(2*TMath::Pi()*(pow(sigma*Amplitude_err,2) + pow(Amplitude*sigma_err,2)));
+ //one_gauss->SetParameter(0,Amplitude);
+ //one_gauss->SetParameter(1,mean);
+ //one_gauss->SetParameter(2,sigma);
+ //one_gauss->SetParError(0,Amplitude_err);
+ //one_gauss->SetParError(1,mean_err);
+ //one_gauss->SetParError(2,sigma_err);
+
+ cout << "************************" << endl;
+ cout << "A= " << A << endl;
+ //cout << "Amplitude= " << Amplitude << "+/-" << Amplitude_err << endl;
+ //cout << "Mean= " << mean << endl;
+ cout << "Sigma= " << sigma << "+/-" << sigma_err << endl;
+
+
+
+ Integral[i] = (12.5/rebin)*Amplitude*sigma*sqrt(2*TMath::Pi());
+ Integral_err[i] = (12.5/rebin)*sqrt(2*TMath::Pi()*(pow(sigma*Amplitude_err,2) + pow(Amplitude*sigma_err,2)));
//Integral_err[i] = sqrt(Integral[i]);
@@ -115,6 +130,7 @@ void Fit(TH1F* hmass,string Energy){
A++;
}
+ double histo_integral = hmass->Integral();
hmass->Write();
total->Write();
@@ -125,9 +141,10 @@ void Fit(TH1F* hmass,string Energy){
// Yield calculation //
A = Amin;
for(int i=0; i<NumberOfA; i++){
- Yield[i] = Integral[i]/total_integral*200;
- gyield->SetPoint(i,A,Yield[i]);
- double yield_err = Integral_err[i]/total_integral*200;
+ //Yield[i] = Integral[i]/total_integral*200;
+ Yield[i] = Integral[i]/histo_integral*200;
+ //double yield_err = Integral_err[i]/total_integral*200;
+ double yield_err = Integral_err[i]/histo_integral*200;
ofile << A << " " << Yield[i] << " " << yield_err << endl;
if(A==99) cout << "A=99 -> Y= " << Yield[i] << endl;
if(A==147) cout << "A=147 -> Y= " << Yield[i] << endl;
@@ -135,7 +152,5 @@ void Fit(TH1F* hmass,string Energy){
}
ofile.close();
- gyield->SetMarkerStyle(8);
- gyield->SetMarkerSize(1);
}
diff --git a/Projects/ana_e850/macro/Mass/GetTotalYield.C b/Projects/ana_e850/macro/Mass/GetTotalYield.C
new file mode 100644
index 000000000..63ce3e31b
--- /dev/null
+++ b/Projects/ana_e850/macro/Mass/GetTotalYield.C
@@ -0,0 +1,143 @@
+void Fit(TH1F* hmass);
+TFile* ofile;
+TGraph* gsigma;
+
+//////////////////////////////////////////////////////
+void GetTotalYield(){
+
+ TFile* ifile = new TFile("histo_mass_all.root","read");
+
+ ofile = new TFile("mass_all_fitted.root","recreate");
+
+ gsigma = new TGraph();
+ TString histo_name = "hmass_all";
+ TH1F* h1 = (TH1F*)ifile->FindObjectAny(histo_name);
+ Fit(h1);
+ ofile->Close();
+
+ TCanvas* c2 = new TCanvas("c2","c2",1200,1200);
+ c2->cd();
+ gsigma->SetMarkerStyle(8);
+ gsigma->Draw("ap");
+
+
+}
+
+//////////////////////////////////////////////////////
+void Fit(TH1F* hmass){
+ int rebin = 1;
+ //hmass->Rebin(rebin);
+ TCanvas* c1 = new TCanvas("c1","c1",1200,1200);
+ c1->cd();
+ hmass->Draw("E1");
+
+ int Amin = 82;
+ int Amax = 153;
+
+ int A = Amin;
+
+ int NumberOfA = Amax - Amin;
+ Double_t para[3*NumberOfA];
+
+ TF1* g[NumberOfA];
+ double Integral[NumberOfA];
+ double Integral_err[NumberOfA];
+ double total_integral = 0;
+ double Yield[NumberOfA];
+
+ for(int i=0; i<NumberOfA; i++){
+ g[i] = new TF1(Form("g%i",i),"gaus",A-0.4,A+0.4);
+ if(i==0){
+ g[i]->SetParameter(0,100);
+ g[i]->SetParameter(1,Amin);
+ g[i]->SetParameter(2,0.3);
+ }
+ else if(i>0){
+ g[i]->SetParameter(0,g[i-1]->GetParameter(0));
+ g[i]->SetParameter(1,A);
+ g[i]->SetParameter(2,g[i-1]->GetParameter(2));
+
+ g[i]->SetParLimits(1,A-0.15,A+0.15);
+ g[i]->SetParLimits(2,0.25,0.40);
+
+ }
+
+ hmass->Fit(g[i],"qBR");
+
+ g[i]->Draw("lsame");
+ g[i]->GetParameters(¶[3*i]);
+
+ A++;
+ }
+
+
+ TString total_func = "gaus(0)";
+ for(int i=1; i<NumberOfA; i++){
+ TString gaus_func = Form("gaus(%i)",3*i);
+ total_func += "+" + gaus_func;
+ }
+
+ TF1* total = new TF1("total",total_func,Amin-0.5, Amin+NumberOfA-0.5);
+ total->SetParameters(para);
+ A = Amin;
+ for(int i=0; i<NumberOfA; i++){
+ total->SetParameter(3*i,g[i]->GetParameter(0));
+ total->SetParameter(3*i+1,g[i]->GetParameter(1));
+ total->SetParameter(3*i+2,g[i]->GetParameter(2));
+
+ total->SetParLimits(3*i+1,A-0.2,A+0.2);
+ total->SetParLimits(3*i+2,g[i]->GetParameter(2)*0.85,g[i]->GetParameter(2)*1.15);
+ A++;
+ }
+ total->SetLineColor(4);
+
+ hmass->Fit(total,"qR");
+
+ A = Amin;
+ for(int i=0; i<NumberOfA; i++){
+
+ double Amplitude = total->GetParameter(3*i);
+ double mean = total->GetParameter(3*i+1);
+ double sigma = total->GetParameter(3*i+2);
+ double Amplitude_err = total->GetParError(3*i);
+ double sigma_err = total->GetParError(3*i+2);
+
+ cout << "A= " << A << endl;
+ cout << "Amplitude= " << Amplitude << " / " << g[i]->GetParameter(0) << endl;
+ cout << "Mean= " << mean << " / " << g[i]->GetParameter(1) << endl;
+ cout << "Sigma= " << sigma << " / " << g[i]->GetParameter(2) << endl;
+
+
+ Integral[i] = (12.5/rebin)*Amplitude*sigma*sqrt(2*TMath::Pi());
+ //Integral_err[i] = sqrt(2*TMath::Pi()*(pow(sigma*Amplitude_err,2) + pow(Amplitude*sigma_err,2)));
+ Integral_err[i] = sqrt(Integral[i]);
+
+ gsigma->SetPoint(i,A,sigma);
+
+ total_integral += Integral[i];
+ A++;
+ }
+
+ double histo_integral = hmass->Integral();
+ hmass->Write();
+ total->Write();
+ gsigma->Write();
+
+ string output_filename = "Total_Yield.dat";
+ ofstream ofile;
+ ofile.open(output_filename.c_str());
+
+ // Yield calculation //
+ A = Amin;
+ for(int i=0; i<NumberOfA; i++){
+ //Yield[i] = Integral[i]/total_integral*200;
+ Yield[i] = Integral[i]/histo_integral*200;
+
+ //double yield_err = Integral_err[i]/total_integral*200;
+ double yield_err = Integral_err[i]/histo_integral*200;
+ ofile << A << " " << Yield[i] << " " << yield_err << endl;
+ A++;
+ }
+ ofile.close();
+
+}
diff --git a/Projects/ana_e850/macro/Mass/GetYieldEvolution.C b/Projects/ana_e850/macro/Mass/GetYieldEvolution.C
new file mode 100644
index 000000000..537c5f7f9
--- /dev/null
+++ b/Projects/ana_e850/macro/Mass/GetYieldEvolution.C
@@ -0,0 +1,179 @@
+void Fit(TH1F* hmass, string Energy);
+
+//////////////////////////////////////////////////////
+void GetYieldEvolution(){
+
+ //TFile* ifile = new TFile("histo_mass_all.root","read");
+ TFile* ifile = new TFile("histo_mass_part2.root","read");
+
+ TH2F* hnorm = new TH2F("hnorm","hnorm",20,0,20,1000,80,160);
+ TH2F* hyield = new TH2F("hyield","hyield",20,0,20,1000,80,160);
+
+ TString histo_name = "hmass_vs_ex";
+ TH2F* h2 = (TH2F*)ifile->FindObjectAny(histo_name);
+
+ h2->Draw("colz");
+
+ double nentries[20];
+ for(int i=0; i<h2->GetNbinsX(); i++){
+ h2->ProjectionY("h2_py",i,i);
+ TH1F* hpy = (TH1F*)gDirectory->FindObjectAny("h2_py");
+
+ nentries[i] = hpy->Integral(1,h2->GetNbinsY());
+
+ cout << i << " " << nentries[i] << endl;
+ }
+
+ // Getting norm
+ for(int i=0; i<hnorm->GetNbinsX(); i++){
+ for(int j=0; j<hnorm->GetNbinsY(); j++){
+ hnorm->SetBinContent(i,j,nentries[i]);
+ }
+ }
+
+ // Filling normalize histogram
+ for(int i=6; i<hnorm->GetNbinsX(); i++){
+ for(int j=0; j<hnorm->GetNbinsY(); j++){
+ double val_orig = h2->GetBinContent(i,j);
+ double val_error = sqrt(val_orig);
+ double val_norm = hnorm->GetBinContent(i,j);
+ double new_val = val_orig/val_norm*200;
+ hyield->SetBinContent(i,j,new_val);
+ hyield->SetBinError(i,j,val_error/val_norm*200);
+ }
+ }
+
+ hyield->Draw("colz");
+
+}
+
+//////////////////////////////////////////////////////
+void Fit(TH1F* hmass,string Energy){
+ int rebin = 2;
+ hmass->Rebin(rebin);
+ hmass->Draw();
+
+ int Amin = 82;
+ int Amax = 153;
+
+ int A = Amin;
+
+ int NumberOfA = Amax - Amin;
+ Double_t para[3*NumberOfA];
+
+ TF1* g[NumberOfA];
+ double Integral[NumberOfA];
+ double Integral_err[NumberOfA];
+ double total_integral = 0;
+ double Yield[NumberOfA];
+
+ for(int i=0; i<NumberOfA; i++){
+ g[i] = new TF1(Form("g%i",i),"gaus",A-0.4,A+0.4);
+ if(i==0){
+ g[i]->SetParameter(0,100);
+ g[i]->SetParameter(1,Amin);
+ g[i]->SetParameter(2,0.3);
+ }
+ else if(i>0){
+ g[i]->SetParameter(0,g[i-1]->GetParameter(0));
+ g[i]->SetParameter(1,g[i-1]->GetParameter(1));
+ g[i]->SetParameter(2,g[i-1]->GetParameter(2));
+
+ g[i]->SetParLimits(1,A-0.15,A+0.15);
+ g[i]->SetParLimits(2,0.1,0.35);
+
+ }
+
+ hmass->Fit(g[i],"qR");
+
+ g[i]->Draw("lsame");
+ g[i]->GetParameters(¶[3*i]);
+
+ A++;
+ }
+
+
+ TString total_func = "gaus(0)";
+ for(int i=1; i<NumberOfA; i++){
+ TString gaus_func = Form("gaus(%i)",3*i);
+ total_func += "+" + gaus_func;
+ }
+
+ TF1* total = new TF1("total",total_func,Amin-0.5, Amin+NumberOfA-0.5);
+
+ A = Amin;
+ for(int i=0; i<NumberOfA; i++){
+ total->SetParameter(3*i,g[i]->GetParameter(0));
+ total->SetParameter(3*i+1,g[i]->GetParameter(1));
+ total->SetParameter(3*i+2,g[i]->GetParameter(2));
+ total->SetParLimits(3*i+1,A-0.15,A+0.15);
+ total->SetParLimits(3*i+2,0.1,0.35);
+ //cout << A << " " << g[i]->GetParameter(2) << endl;
+ A++;
+ }
+ total->SetParameters(para);
+ total->SetLineColor(4);
+
+ hmass->Fit(total,"Rq");
+
+ A = Amin;
+ TF1* one_gauss = new TF1("one_gaus","gaus");
+ for(int i=0; i<NumberOfA; i++){
+
+ double Amplitude = total->GetParameter(3*i);
+ double mean = total->GetParameter(3*i+1);
+ double sigma = total->GetParameter(3*i+2);
+ double Amplitude_err = total->GetParError(3*i);
+ double mean_err = total->GetParError(3*i+1);
+ double sigma_err = total->GetParError(3*i+2);
+
+ one_gauss->SetParameter(0,Amplitude);
+ one_gauss->SetParameter(1,mean);
+ one_gauss->SetParameter(2,sigma);
+ one_gauss->SetParError(0,Amplitude_err);
+ one_gauss->SetParError(1,mean_err);
+ one_gauss->SetParError(2,sigma_err);
+
+ //cout << "A= " << A << endl;
+ //cout << "Amplitude= " << Amplitude << "+/-" << Amplitude_err << endl;
+ //cout << "Mean= " << mean << endl;
+ //cout << "Sigma= " << sigma << "+/-" << sigma_err << endl;
+
+
+
+ Integral[i] = (12.5/rebin)*Amplitude*sigma*sqrt(2*TMath::Pi());
+ Integral_err[i] = (12.5/rebin)*sqrt(2*TMath::Pi()*(pow(sigma*Amplitude_err,2) + pow(Amplitude*sigma_err,2)));
+ //Integral_err[i] = sqrt(Integral[i]);
+
+ cout << "**************** Compare integral " << endl;
+ cout << Integral[i] << " / " << (12.5/rebin)*one_gauss->Integral(A-1,A+1) << endl;
+ cout << Integral_err[i] << " / " << one_gauss->IntegralError(A-1,A+1) << endl;
+
+ total_integral += Integral[i];
+ A++;
+ }
+
+ double histo_integral = hmass->Integral();
+ hmass->Write();
+ total->Write();
+
+ string output_filename = "Mass_Yield_" + Energy + ".dat";
+ ofstream ofile;
+ ofile.open(output_filename.c_str());
+
+ // Yield calculation //
+ A = Amin;
+ for(int i=0; i<NumberOfA; i++){
+ Yield[i] = Integral[i]/total_integral*200;
+ //Yield[i] = Integral[i]/histo_integral*200;
+ double yield_err = Integral_err[i]/total_integral*200;
+ //double yield_err = Integral_err[i]/histo_integral*200;
+ ofile << A << " " << Yield[i] << " " << yield_err << endl;
+ if(A==99) cout << "A=99 -> Y= " << Yield[i] << endl;
+ if(A==147) cout << "A=147 -> Y= " << Yield[i] << endl;
+ A++;
+ }
+ ofile.close();
+
+
+}
diff --git a/Projects/ana_e850/macro/Mass/PlotYieldEvolution.C b/Projects/ana_e850/macro/Mass/PlotYieldEvolution.C
index eac2ab885..41c25a808 100644
--- a/Projects/ana_e850/macro/Mass/PlotYieldEvolution.C
+++ b/Projects/ana_e850/macro/Mass/PlotYieldEvolution.C
@@ -1,6 +1,7 @@
-void PlotYieldEvolution(int A_asked=147){
- int Energy[10] = {6,7,8,9,10,11,12,13,14,15};
+void PlotYieldEvolution(int A_asked=99){
+ int Energy[11] = {6,7,8,9,10,11,12,13,14,15,16};
TGraphErrors* gevol = new TGraphErrors();
+ TGraphErrors* g_relative = new TGraphErrors();
for(int i = 0; i<10; i++){
ifstream ifile;
@@ -10,15 +11,26 @@ void PlotYieldEvolution(int A_asked=147){
int A;
double yield;
double yield_err;
+ double Sn = 6.534;
while(ifile>>A>>yield>>yield_err){
if(A==A_asked){
- gevol->SetPoint(i,Energy[i]-5.5,yield);
+ gevol->SetPoint(i,Energy[i]-Sn,yield);
gevol->SetPointError(i,0,yield_err);
+
+ g_relative->SetPoint(i,Energy[i]-5.5,yield_err/yield);
}
}
}
+ TH2F* h2 = new TH2F("h2","h2",100,-1,10,100,1,8);
+ TCanvas* c1 = new TCanvas("c1","c1",800,1600);
+ c1->Divide(1,2);
+ c1->cd(1);
gevol->SetMarkerStyle(8);
gevol->SetMarkerSize(1);
- gevol->Draw("ap");
+ h2->Draw();
+ gevol->Draw("psame");
+ c1->cd(2);
+ g_relative->SetMarkerStyle(8);
+ g_relative->Draw("ap");
}
diff --git a/Projects/ana_e850/macro/Mass/check_gamma.C b/Projects/ana_e850/macro/Mass/check_gamma.C
index aa2a61426..1c1c09590 100644
--- a/Projects/ana_e850/macro/Mass/check_gamma.C
+++ b/Projects/ana_e850/macro/Mass/check_gamma.C
@@ -167,5 +167,7 @@ void check_gamma(string part="part1")
hg100->Write();
hg98->Write();
+ hg10Be->Write();
+
ofile->Close();
}
--
GitLab