From 02750723f7ecf3f8e7cb9afa99530714b2e85d87 Mon Sep 17 00:00:00 2001
From: Morfouace <pierre.morfouace@gmail.com>
Date: Mon, 2 Oct 2023 15:06:26 +0200
Subject: [PATCH] Updating e850 macro
---
Projects/e850/macro/Mass/FillHistoMass.C | 139 ++++++++++++++++++
Projects/e850/macro/Mass/FitMass.C | 123 ++++++++++++++++
Projects/e850/macro/Mass/PlotYieldEvolution.C | 24 +++
Projects/e850/macro/Mass/cut/Draw_cut.C | 35 +++++
4 files changed, 321 insertions(+)
create mode 100644 Projects/e850/macro/Mass/FillHistoMass.C
create mode 100644 Projects/e850/macro/Mass/FitMass.C
create mode 100644 Projects/e850/macro/Mass/PlotYieldEvolution.C
create mode 100644 Projects/e850/macro/Mass/cut/Draw_cut.C
diff --git a/Projects/e850/macro/Mass/FillHistoMass.C b/Projects/e850/macro/Mass/FillHistoMass.C
new file mode 100644
index 000000000..f2d303d66
--- /dev/null
+++ b/Projects/e850/macro/Mass/FillHistoMass.C
@@ -0,0 +1,139 @@
+TChain* chain;
+TCutG* cut10Be[8];
+
+/////////////////////////////////////////////////////////////////
+void LoadCut(string part){
+ string path = "cut/"+ part + "/cut10Be.root";
+ TFile* fcut = new TFile(path.c_str());
+ for(int i=0; i<8; i++){
+ cut10Be[i] = (TCutG*)fcut->Get(Form("cut10Be_det%i",i+1));
+ }
+}
+
+/////////////////////////////////////////////////////////////////
+void FillHistoMass(string part="part2"){
+
+ LoadCut(part);
+
+ chain = new TChain("PhysicsTree");
+ if(part=="part1"){
+ chain->Add("../../root/analysis/run_28.root");
+ chain->Add("../../root/analysis/run_29.root");
+ chain->Add("../../root/analysis/run_30.root");
+ chain->Add("../../root/analysis/run_31.root");
+ chain->Add("../../root/analysis/run_32.root");
+ chain->Add("../../root/analysis/run_36.root");
+ chain->Add("../../root/analysis/run_37.root");
+ chain->Add("../../root/analysis/run_39.root");
+ chain->Add("../../root/analysis/run_40.root");
+ chain->Add("../../root/analysis/run_41.root");
+ chain->Add("../../root/analysis/run_42.root");
+ chain->Add("../../root/analysis/run_43.root");
+ chain->Add("../../root/analysis/run_44.root");
+ chain->Add("../../root/analysis/run_45.root");
+ chain->Add("../../root/analysis/run_46.root");
+ chain->Add("../../root/analysis/run_48.root");
+ chain->Add("../../root/analysis/run_49.root");
+ chain->Add("../../root/analysis/run_50.root");
+ chain->Add("../../root/analysis/run_51.root");
+ chain->Add("../../root/analysis/run_52.root");
+ chain->Add("../../root/analysis/run_53.root");
+ }
+ else if(part=="part2"){
+ chain->Add("../../root/analysis/run_55.root");
+ chain->Add("../../root/analysis/run_56.root");
+ chain->Add("../../root/analysis/run_57.root");
+ chain->Add("../../root/analysis/run_58.root");
+ chain->Add("../../root/analysis/run_59.root");
+ chain->Add("../../root/analysis/run_60.root");
+ chain->Add("../../root/analysis/run_61.root");
+ chain->Add("../../root/analysis/run_62.root");
+ chain->Add("../../root/analysis/run_63.root");
+ chain->Add("../../root/analysis/run_64.root");
+ chain->Add("../../root/analysis/run_65.root");
+ chain->Add("../../root/analysis/run_66.root");
+ chain->Add("../../root/analysis/run_67.root");
+ chain->Add("../../root/analysis/run_68.root");
+ chain->Add("../../root/analysis/run_69.root");
+ }
+ double DeltaEcorr;
+ double Elab;
+ double Mass;
+ double Ex240Pu;
+ double fTS_TMW;
+ int Telescope;
+
+ chain->SetBranchStatus("DeltaEcorr","true");
+ chain->SetBranchAddress("DeltaEcorr",&DeltaEcorr);
+ chain->SetBranchStatus("Elab","true");
+ chain->SetBranchAddress("Elab",&Elab);
+ chain->SetBranchStatus("Mass","true");
+ chain->SetBranchAddress("Mass",&Mass);
+ chain->SetBranchStatus("Ex240Pu","true");
+ chain->SetBranchAddress("Ex240Pu",&Ex240Pu);
+ chain->SetBranchStatus("fTS_TMW","true");
+ chain->SetBranchAddress("fTS_TMW",&fTS_TMW);
+ chain->SetBranchStatus("Telescope","true");
+ chain->SetBranchAddress("Telescope",&Telescope);
+
+
+
+ TH1F *hmass[11];
+ TH2F *hPID[8];
+ for(int i=0; i<11; i++){
+ TString histo_name = Form("hmass_%iMeV",i+6);
+ hmass[i] = new TH1F(histo_name, histo_name,1000,80,160);
+ }
+ for(int i=0; i<8; i++){
+ TString histo_name = Form("PID_det%i",i+1);
+ hPID[i] = new TH2F(histo_name,histo_name,1000,0,200,1000,0,60);
+ }
+ //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);
+
+
+ int nentries = chain->GetEntries();
+ cout << "**** Number of entries to be treated: " << nentries << " ****" << endl;
+ for(int i=0; i<nentries; i++){
+ chain->GetEntry(i);
+ if(i%1000000==0){
+ cout << "\033[34m\r Processing tree..." << (double)i/nentries*100 << "\% done" << flush;
+ }
+
+ if(fTS_TMW>0){
+
+ if(Telescope>0) hPID[Telescope-1]->Fill(Elab,DeltaEcorr);
+
+ if(cut10Be[Telescope-1]->IsInside(Elab,DeltaEcorr)){
+ hEx->Fill(Ex240Pu);
+ hMassAll->Fill(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);
+ else if(Ex240Pu>8.5 && Ex240Pu<9.5) hmass[3]->Fill(Mass);
+ else if(Ex240Pu>9.5 && Ex240Pu<10.5) hmass[4]->Fill(Mass);
+ else if(Ex240Pu>10.5 && Ex240Pu<11.5) hmass[5]->Fill(Mass);
+ else if(Ex240Pu>11.5 && Ex240Pu<12.5) hmass[6]->Fill(Mass);
+ else if(Ex240Pu>12.5 && Ex240Pu<13.5) hmass[7]->Fill(Mass);
+ else if(Ex240Pu>13.5 && Ex240Pu<14.5) hmass[8]->Fill(Mass);
+ else if(Ex240Pu>14.5 && Ex240Pu<15.5) hmass[9]->Fill(Mass);
+ else if(Ex240Pu>15.5 && Ex240Pu<16.5) hmass[10]->Fill(Mass);
+ }
+ }
+ }
+
+ TString filename = Form("histo_mass_%s.root",part.c_str());
+ TFile *ofile = new TFile(filename,"recreate");
+ for(int i=0; i<11; i++){
+ hmass[i]->Write();
+ }
+ hEx->Write();
+ hMassAll->Write();
+ for(int i=0; i<8; i++)
+ hPID[i]->Write();
+ ofile->Close();
+
+
+
+}
diff --git a/Projects/e850/macro/Mass/FitMass.C b/Projects/e850/macro/Mass/FitMass.C
new file mode 100644
index 000000000..f4f4b2978
--- /dev/null
+++ b/Projects/e850/macro/Mass/FitMass.C
@@ -0,0 +1,123 @@
+void Fit(TH1F* hmass, string Energy);
+TFile* ofile;
+
+//////////////////////////////////////////////////////
+void FitMass(){
+
+ TFile* ifile = new TFile("histo_mass_all.root","read");
+
+ ofile = new TFile("histo_mass_fitted.root","recreate");
+
+ for(int i=0; i<10; i++){
+ int Ex = i+6;
+ string Energy = to_string(Ex) + "MeV";
+ TString histo_name = Form("hmass_%iMeV",i+6);
+ TH1F* h1 = (TH1F*)ifile->FindObjectAny(histo_name);
+ Fit(h1,Energy);
+ }
+
+ ofile->Close();
+}
+
+//////////////////////////////////////////////////////
+void Fit(TH1F* hmass,string Energy){
+ hmass->Rebin(2);
+ hmass->Draw();
+
+ TGraph* gyield = new TGraph();
+ 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);
+ g[i]->SetParLimits(0,0.1,20000);
+ g[i]->SetParLimits(1,A-0.1,A+0.1);
+ g[i]->SetParLimits(2,0.1,0.45);
+
+ 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->SetParLimits(3*i+1,A-0.05,A+0.05);
+ total->SetParLimits(3*i+2,0.25,0.35);
+ A++;
+ }
+ total->SetParameters(para);
+ total->SetLineColor(4);
+
+ hmass->Fit(total,"Rq");
+
+ 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 << 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)));
+ //Integral_err[i] = sqrt(Integral[i]);
+
+
+ total_integral += Integral[i];
+ A++;
+ }
+
+ 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;
+ gyield->SetPoint(i,A,Yield[i]);
+ double yield_err = Integral_err[i]/total_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();
+
+ gyield->SetMarkerStyle(8);
+ gyield->SetMarkerSize(1);
+
+}
diff --git a/Projects/e850/macro/Mass/PlotYieldEvolution.C b/Projects/e850/macro/Mass/PlotYieldEvolution.C
new file mode 100644
index 000000000..eac2ab885
--- /dev/null
+++ b/Projects/e850/macro/Mass/PlotYieldEvolution.C
@@ -0,0 +1,24 @@
+void PlotYieldEvolution(int A_asked=147){
+ int Energy[10] = {6,7,8,9,10,11,12,13,14,15};
+ TGraphErrors* gevol = new TGraphErrors();
+
+ for(int i = 0; i<10; i++){
+ ifstream ifile;
+ string input_filename = "Mass_Yield_" + to_string(Energy[i]) + "MeV.dat";
+ ifile.open(input_filename.c_str());
+
+ int A;
+ double yield;
+ double yield_err;
+ while(ifile>>A>>yield>>yield_err){
+ if(A==A_asked){
+ gevol->SetPoint(i,Energy[i]-5.5,yield);
+ gevol->SetPointError(i,0,yield_err);
+ }
+ }
+ }
+
+ gevol->SetMarkerStyle(8);
+ gevol->SetMarkerSize(1);
+ gevol->Draw("ap");
+}
diff --git a/Projects/e850/macro/Mass/cut/Draw_cut.C b/Projects/e850/macro/Mass/cut/Draw_cut.C
new file mode 100644
index 000000000..cf336bf20
--- /dev/null
+++ b/Projects/e850/macro/Mass/cut/Draw_cut.C
@@ -0,0 +1,35 @@
+
+//////////////////////////////////////////
+void Draw_cut()
+{
+
+ TChain* chain = new TChain("PhysicsTree");
+ chain->Add("../../../root/analysis/run_55.root");
+ //TFile* histofile = new TFile("histo/histo_EQ1_lg.root");
+ TFile* cutfile = new TFile("cut10Be.root","update");
+
+ TCanvas* c1 = new TCanvas("c1","c1",1200,1200);
+
+
+ for(int det_number=1; det_number<9; det_number++){
+ TString to_draw = Form("DeltaEcorr:Elab>>h%i(1000,0,200,1000,0,60)",det_number);
+ TString condition = Form("fTS_TMW>0 && Telescope==%i",det_number);
+
+ chain->Draw(to_draw,condition,"colz",5e7);
+
+ TH2F *h = (TH2F*)gDirectory->FindObjectAny(Form("h%i",det_number));
+ h->Draw("colz");
+
+ TCutG* cutg;
+ cutg = (TCutG*)gPad->WaitPrimitive("CUTG","CutG");
+ cutg->SetName(Form("cut10Be_det%i",det_number));
+
+ cutfile->cd();
+ cutg->Write();
+ }
+ //histofile->Close();
+ cutfile->Close();
+
+}
+
+/////////////////////////////////////////
--
GitLab