diff --git a/NPLib/scripts/build_dict.sh.in b/NPLib/scripts/build_dict.sh.in
index e032afeff9cffb70e6746e129e76136542b319e4..7e810fcbb83cf4265db592ef3f901ef0f31561c8 100755
--- a/NPLib/scripts/build_dict.sh.in
+++ b/NPLib/scripts/build_dict.sh.in
@@ -39,6 +39,12 @@ fi
# Get the Major version
version_major_string="${version%.*}"
+# Since 6.30 they changed the numbering scheme from 6.XX/XX
+# to 6.XX.XX
+num_char=${#version_major_string}
+if [ $num_char -gt 1 ]
+then version_major_string="${version%.*.*}"
+fi
version_major=$(($version_major_string))
# if before version 4, exit
diff --git a/Projects/e870/Analysis.cxx b/Projects/e870/Analysis.cxx
index af55fe7ba7ff505fa1267fe2f1d7178ab39cf748..a3345f779b1543846b1114d8a7ffbaf9cc231b40 100644
--- a/Projects/e870/Analysis.cxx
+++ b/Projects/e870/Analysis.cxx
@@ -21,7 +21,7 @@
*****************************************************************************/
#include <iostream>
using namespace std;
-
+#include "NPVAnalysis.h"
#include "Analysis.h"
#include "NPAnalysisFactory.h"
#include "NPDetectorManager.h"
@@ -91,9 +91,6 @@ void Analysis::Init() {
CsI_E_M2 = 0;
Energy = 0;
ThetaGDSurface = 0;
- X = 0;
- Y = 0;
- Z = 0;
dE = 0;
BeamDirection = TVector3(0, 0, 1);
nbTrack = 0;
@@ -128,19 +125,24 @@ void Analysis::TreatEvent() {
int TelescopeNumber = M2->TelescopeNumber[countMust2];
/************************************************/
+
// Part 1 : Impact Angle
+ //THIS IS THE NEW PART
+
ThetaM2Surface = 0;
ThetaNormalTarget = 0;
TVector3 HitDirection = M2->GetPositionOfInteraction(countMust2) - BeamImpact;
+ TVector3 Coords = M2->GetPositionOfInteraction(countMust2);
double Theta = HitDirection.Angle(BeamDirection);
-
- X = M2->GetPositionOfInteraction(countMust2).X();
- Y = M2->GetPositionOfInteraction(countMust2).Y();
- Z = M2->GetPositionOfInteraction(countMust2).Z();
+ double Phi = HitDirection.Phi();
+ X.push_back(Coords.x());
+ Y.push_back(Coords.y());
+ Z.push_back(Coords.z());
ThetaM2Surface = HitDirection.Angle(-M2->GetTelescopeNormal(countMust2));
ThetaNormalTarget = HitDirection.Angle(TVector3(0, 0, 1));
+
/************************************************/
/************************************************/
@@ -167,9 +169,12 @@ void Analysis::TreatEvent() {
Energy = LightTarget.EvaluateInitialEnergy(Energy, TargetThickness * 0.5, Theta);
// What is written in the tree
+
ThetaLab.push_back(Theta / deg);
+ PhiLab.push_back(Phi / deg);
Ex.push_back(reaction.ReconstructRelativistic(Energy, Theta));
ELab.push_back(Energy);
+
/************************************************/
} // end loop MUST2
@@ -183,11 +188,15 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("Ex", &Ex);
RootOutput::getInstance()->GetTree()->Branch("ELab", &ELab);
RootOutput::getInstance()->GetTree()->Branch("ThetaLab", &ThetaLab);
+ RootOutput::getInstance()->GetTree()->Branch("PhiLab", &PhiLab);
RootOutput::getInstance()->GetTree()->Branch("ThetaCM", &ThetaCM, "ThetaCM/D");
RootOutput::getInstance()->GetTree()->Branch("Run", &Run, "Run/I");
- RootOutput::getInstance()->GetTree()->Branch("X", &X, "X/D");
- RootOutput::getInstance()->GetTree()->Branch("Y", &Y, "Y/D");
- RootOutput::getInstance()->GetTree()->Branch("Z", &Z, "Z/D");
+ // RootOutput::getInstance()->GetTree()->Branch("X", &X, "X/D"); Like before
+ // RootOutput::getInstance()->GetTree()->Branch("Y", &Y, "Y/D");
+ // RootOutput::getInstance()->GetTree()->Branch("Z", &Z, "Z/D");
+ RootOutput::getInstance()->GetTree()->Branch("X", &X);
+ RootOutput::getInstance()->GetTree()->Branch("Y", &Y);
+ RootOutput::getInstance()->GetTree()->Branch("Z", &Z);
RootOutput::getInstance()->GetTree()->Branch("dE", &dE, "dE/D");
if (!simulation) {
}
@@ -222,12 +231,16 @@ void Analysis::ReInitValue() {
EDC = -1000;
BeamEnergy = -1000;
ThetaCM = -1000;
- X = -1000;
- Y = -1000;
- Z = -1000;
+ // X = -1000;
+ // Y = -1000;
+ // Z = -1000;
+ X.clear();
+ Y.clear();
+ Z.clear();
dE = -1000;
ELab.clear();
ThetaLab.clear();
+ PhiLab.clear();
Ex.clear();
}
diff --git a/Projects/e870/Analysis.h b/Projects/e870/Analysis.h
index 753d9fcdbfb81f25e6efff8ba3727f67561884b3..afe2d339f85511b3826f370f188c4b3248655e80 100644
--- a/Projects/e870/Analysis.h
+++ b/Projects/e870/Analysis.h
@@ -58,6 +58,7 @@ class Analysis: public NPL::VAnalysis{
double EDC;
std::vector<double> ELab;
std::vector<double> ThetaLab;
+ std::vector<double> PhiLab;
double ThetaCM;
double OriginalELab;
double OriginalThetaLab;
@@ -93,9 +94,9 @@ class Analysis: public NPL::VAnalysis{
double BeamEnergy;
double ThetaGDSurface ;
- double X ;
- double Y ;
- double Z ;
+ std::vector<double> X ;
+ std::vector<double> Y ;
+ std::vector<double> Z ;
// Vamos Branches
unsigned long long int LTS;
diff --git a/Projects/e870/DetectorConfiguration/MUST2_E805.detector b/Projects/e870/DetectorConfiguration/MUST2_E805.detector
index 4c78f7563ad7c4d4783c632bee6c0a61c96cc24f..bc14fd00cb6f8c2b1fd2902b90aac644f6608fd6 100644
--- a/Projects/e870/DetectorConfiguration/MUST2_E805.detector
+++ b/Projects/e870/DetectorConfiguration/MUST2_E805.detector
@@ -58,3 +58,18 @@ M2Telescope
VIS= all
CsIOffset= 1
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+CATSDetector
+ CATSNumber = 1
+ X1_Y1= 43.71 -41.76 -1587.1 mm
+ X28_Y1= -27.41 -41.76 -1587.1 mm
+ X1_Y28= 43.71 29.36 -1587.1 mm
+ X28_Y28= -27.41 29.36 -1587.1 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+CATSDetector
+ CATSNumber = 2
+ X1_Y1= 33.94 -37.06 -1090.1 mm
+ X28_Y1= -37.18 -37.06 -1090.1 mm
+ X1_Y28= 33.94 34.06 -1090.1 mm
+ X28_Y28= -37.18 34.06 -1090.1 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
diff --git a/Projects/e870/Makefile b/Projects/e870/Makefile
index 719c7c7a4e2b71e65030d87930757eb33ae68a19..223d5e7e565eed065829403cf49c6917558f3848 100644
--- a/Projects/e870/Makefile
+++ b/Projects/e870/Makefile
@@ -1,44 +1,35 @@
# CMAKE generated file: DO NOT EDIT!
-# Generated by "Unix Makefiles" Generator, CMake Version 3.26
+# Generated by "Unix Makefiles" Generator, CMake Version 3.16
# Default target executed when no arguments are given to make.
default_target: all
+
.PHONY : default_target
# Allow only one "make -f Makefile2" at a time, but pass parallelism.
.NOTPARALLEL:
+
#=============================================================================
# Special targets provided by cmake.
# Disable implicit rules so canonical targets will work.
.SUFFIXES:
-# Disable VCS-based implicit rules.
-% : %,v
-
-# Disable VCS-based implicit rules.
-% : RCS/%
-
-# Disable VCS-based implicit rules.
-% : RCS/%,v
-
-# Disable VCS-based implicit rules.
-% : SCCS/s.%
-# Disable VCS-based implicit rules.
-% : s.%
+# Remove some rules from gmake that .SUFFIXES does not remove.
+SUFFIXES =
.SUFFIXES: .hpux_make_needs_suffix_list
-# Command-line flag to silence nested $(MAKE).
-$(VERBOSE)MAKESILENT = -s
-#Suppress display of executed commands.
+# Suppress display of executed commands.
$(VERBOSE).SILENT:
+
# A target that is always out of date.
cmake_force:
+
.PHONY : cmake_force
#=============================================================================
@@ -48,67 +39,70 @@ cmake_force:
SHELL = /bin/sh
# The CMake executable.
-CMAKE_COMMAND = /usr/local/Cellar/cmake/3.26.4/bin/cmake
+CMAKE_COMMAND = /usr/bin/cmake
# The command to remove a file.
-RM = /usr/local/Cellar/cmake/3.26.4/bin/cmake -E rm -f
+RM = /usr/bin/cmake -E remove -f
# Escaping for special characters.
EQUALS = =
# The top-level source directory on which CMake was run.
-CMAKE_SOURCE_DIR = /Users/valerian/Software/nptool_gitlab/nptool/Projects/e870
+CMAKE_SOURCE_DIR = /home/tania/Software/nptool/Projects/e870
# The top-level build directory on which CMake was run.
-CMAKE_BINARY_DIR = /Users/valerian/Software/nptool_gitlab/nptool/Projects/e870
+CMAKE_BINARY_DIR = /home/tania/Software/nptool/Projects/e870
#=============================================================================
# Targets provided globally by CMake.
-# Special rule for the target edit_cache
-edit_cache:
- @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake cache editor..."
- /usr/local/Cellar/cmake/3.26.4/bin/ccmake -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
-.PHONY : edit_cache
-
-# Special rule for the target edit_cache
-edit_cache/fast: edit_cache
-.PHONY : edit_cache/fast
-
# Special rule for the target rebuild_cache
rebuild_cache:
@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake to regenerate build system..."
- /usr/local/Cellar/cmake/3.26.4/bin/cmake --regenerate-during-build -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
+ /usr/bin/cmake -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
.PHONY : rebuild_cache
# Special rule for the target rebuild_cache
rebuild_cache/fast: rebuild_cache
+
.PHONY : rebuild_cache/fast
+# Special rule for the target edit_cache
+edit_cache:
+ @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake cache editor..."
+ /usr/bin/ccmake -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
+.PHONY : edit_cache
+
+# Special rule for the target edit_cache
+edit_cache/fast: edit_cache
+
+.PHONY : edit_cache/fast
+
# The main all target
all: cmake_check_build_system
- $(CMAKE_COMMAND) -E cmake_progress_start /Users/valerian/Software/nptool_gitlab/nptool/Projects/e870/CMakeFiles /Users/valerian/Software/nptool_gitlab/nptool/Projects/e870//CMakeFiles/progress.marks
- $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 all
- $(CMAKE_COMMAND) -E cmake_progress_start /Users/valerian/Software/nptool_gitlab/nptool/Projects/e870/CMakeFiles 0
+ $(CMAKE_COMMAND) -E cmake_progress_start /home/tania/Software/nptool/Projects/e870/CMakeFiles /home/tania/Software/nptool/Projects/e870/CMakeFiles/progress.marks
+ $(MAKE) -f CMakeFiles/Makefile2 all
+ $(CMAKE_COMMAND) -E cmake_progress_start /home/tania/Software/nptool/Projects/e870/CMakeFiles 0
.PHONY : all
# The main clean target
clean:
- $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 clean
+ $(MAKE) -f CMakeFiles/Makefile2 clean
.PHONY : clean
# The main clean target
clean/fast: clean
+
.PHONY : clean/fast
# Prepare targets for installation.
preinstall: all
- $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 preinstall
+ $(MAKE) -f CMakeFiles/Makefile2 preinstall
.PHONY : preinstall
# Prepare targets for installation.
preinstall/fast:
- $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 preinstall
+ $(MAKE) -f CMakeFiles/Makefile2 preinstall
.PHONY : preinstall/fast
# clear depends
@@ -121,36 +115,39 @@ depend:
# Build rule for target.
NPAnalysis: cmake_check_build_system
- $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 NPAnalysis
+ $(MAKE) -f CMakeFiles/Makefile2 NPAnalysis
.PHONY : NPAnalysis
# fast build rule for target.
NPAnalysis/fast:
- $(MAKE) $(MAKESILENT) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/build
+ $(MAKE) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/build
.PHONY : NPAnalysis/fast
Analysis.o: Analysis.cxx.o
+
.PHONY : Analysis.o
# target to build an object file
Analysis.cxx.o:
- $(MAKE) $(MAKESILENT) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/Analysis.cxx.o
+ $(MAKE) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/Analysis.cxx.o
.PHONY : Analysis.cxx.o
Analysis.i: Analysis.cxx.i
+
.PHONY : Analysis.i
# target to preprocess a source file
Analysis.cxx.i:
- $(MAKE) $(MAKESILENT) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/Analysis.cxx.i
+ $(MAKE) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/Analysis.cxx.i
.PHONY : Analysis.cxx.i
Analysis.s: Analysis.cxx.s
+
.PHONY : Analysis.s
# target to generate assembly for a file
Analysis.cxx.s:
- $(MAKE) $(MAKESILENT) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/Analysis.cxx.s
+ $(MAKE) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/Analysis.cxx.s
.PHONY : Analysis.cxx.s
# Help Target
@@ -159,8 +156,8 @@ help:
@echo "... all (the default if no target is provided)"
@echo "... clean"
@echo "... depend"
- @echo "... edit_cache"
@echo "... rebuild_cache"
+ @echo "... edit_cache"
@echo "... NPAnalysis"
@echo "... Analysis.o"
@echo "... Analysis.i"
diff --git a/Projects/e870/macros/DrawExInvariant.cxx b/Projects/e870/macros/DrawExInvariant.cxx
index 0d191812285a13367b9dcfdda01464145e69b49d..c624a907a909411fbe7cbb646db2af37e78d5106 100644
--- a/Projects/e870/macros/DrawExInvariant.cxx
+++ b/Projects/e870/macros/DrawExInvariant.cxx
@@ -1,166 +1,187 @@
-void DrawDalitz() {
+void DrawExInvariant() {
////////////////////////////////////////////////////////////////////////////////
- TLorentzVector LV_alpha, LV_proton, LV_fragment;
- double mu = 931.5; // MeV/c2
- double m_proton = mu + 7.289; // Mev/c2
- double m_alpha = 4 * mu + 2.425; // Mev/c2
- double m_fragment = 44 * mu - 3.755; // Mev/c2
- double m_5Li = 5 * mu + 11.68; // MeV/c2
+ TLorentzVector LV_alpha, LV_7Li, LV_p;
+ double mu = 931.5; // MeV/c2
+ double m_alpha = 4 * mu + 2.425; // Mev/c2
+ double m_7Li = 7 * mu + 14.91; // MeV/c2
+ double m_10Be = 10 * mu + 12.61; // MeV/c2
+ double mp = mu + 7.289; // MeV/c2
+ // double m_5Li = 5 * mu + 11.68; // MeV/c2
- TChain* t = new TChain("PhysicsTree"); // SimulatedTree
- t->Add("NPRootA/NPr0371_*a.root"); // Name of simulated file
+ TChain* t = new TChain("PhysicsTree");
+ t->Add("../../Outputs/Analysis/Ana10Bepalpha_38AMeV_CH2.root"); // Name of simulated file
// YOUR ALPHA AND Li CUT
- TFile* fcuta = new TFile("CUT_e805/CUT_alpha_must.root", "read");
- TCutG* cuta = (TCutG*)gDirectory->FindObjectAny("CUT_alpha_must");
- TFile* fcutp = new TFile("CUT_e805/CUT_proton_must.root", "read");
- TCutG* cutp = (TCutG*)gDirectory->FindObjectAny("CUT_proton_must");
+ TFile* fcuta = new TFile("cuts/CUT_alpha.root", "read");
+ TCutG* cuta = (TCutG*)gDirectory->FindObjectAny("CUT_alpha");
+ TFile* fcutli = new TFile("cuts/CUT_7Li.root", "read");
+ TCutG* cutli = (TCutG*)gDirectory->FindObjectAny("CUT_7Li");
////////////////////////////////////////////////////////////////////////////////
TMust2Physics* M2 = new TMust2Physics();
std::vector<double>* ELab = 0;
std::vector<double>* ThetaLab = 0;
- std::vector<double>* M2_X = 0;
- std::vector<double>* M2_Y = 0;
- std::vector<double>* M2_Z = 0;
+ std::vector<double>* PhiLab = 0;
+ std::vector<double>* Ex = 0;
TBranch* b_ELab;
TBranch* b_ThetaLab;
- TBranch* b_M2_X;
- TBranch* b_M2_Y;
- TBranch* b_M2_Z;
-
- float CATS1_X, CATS1_Y, CATS2_X, CATS2_Y, Xf, Yf;
+ TBranch* b_PhiLab;
+ TBranch* b_X;
+ TBranch* b_Y;
+ TBranch* b_Z;
+ TBranch* b_Ex;
t->SetBranchStatus("*", false);
t->SetBranchStatus("MUST2", true);
t->SetBranchAddress("MUST2", &M2);
- t->SetBranchStatus("M2_ELab", true);
- t->SetBranchAddress("M2_ELab", &ELab, &b_ELab);
- t->SetBranchStatus("M2_ThetaLab", true);
- t->SetBranchAddress("M2_ThetaLab", &ThetaLab, &b_ThetaLab);
- t->SetBranchStatus("CATS1_X", "true");
- t->SetBranchAddress("CATS1_X", &CATS1_X);
- t->SetBranchStatus("CATS1_Y", "true");
- t->SetBranchAddress("CATS1_Y", &CATS1_Y);
- t->SetBranchStatus("CATS2_X", "true");
- t->SetBranchAddress("CATS2_X", &CATS2_X);
- t->SetBranchStatus("CATS2_Y", "true");
- t->SetBranchAddress("CATS2_Y", &CATS2_Y);
- t->SetBranchStatus("M2_X", "true");
- t->SetBranchAddress("M2_X", &M2_X, &b_M2_X);
- t->SetBranchStatus("M2_Y", "true");
- t->SetBranchAddress("M2_Y", &M2_Y, &b_M2_Y);
- t->SetBranchStatus("M2_Z", "true");
- t->SetBranchAddress("M2_Z", &M2_Z, &b_M2_Z);
-
- TH1F* hExTot = new TH1F("hExTot", "hExTot", 1000, -100, 100);
- TH1F* hEx5Li = new TH1F("hEx5Li", "hEx5Li", 1000, -100, 100);
-
+ t->SetBranchStatus("ELab", true);
+ t->SetBranchAddress("ELab", &ELab, &b_ELab);
+ t->SetBranchStatus("ThetaLab", true);
+ t->SetBranchAddress("ThetaLab", &ThetaLab, &b_ThetaLab);
+ t->SetBranchStatus("PhiLab", true);
+ t->SetBranchAddress("PhiLab", &PhiLab, &b_PhiLab);
+ t->SetBranchStatus("Ex", "true");
+ t->SetBranchAddress("Ex", &Ex, &b_Ex);
+
+ TH1F* hExInv = new TH1F("hExInv", "hExInv", 1000, -20, 20);
+ TH1F* hExMM = new TH1F("hExMM", "hExMM", 1000, -20, 20);
TH2F* hE1E2 = new TH2F("hE1E2", "hE1E2", 1000, 0, 500, 1000, 0, 500);
+ TH2F* hELabThetaLab = new TH2F("hELabThetLab", "hELabThetaLab", 1000, 0, 50, 1000, 0, 500);
////////////////////////////////////////////////////////////////////////////////
- // int n_entries = t->GetEntries();
- int n_entries = 1e7;
+ int n_entries = t->GetEntries();
+ // int n_entries = 1e7;
for (unsigned int i = 0; i < n_entries; i++) {
t->GetEntry(i);
- double e_alpha;
- double theta_alpha;
- double phi_alpha;
-
- double e_proton;
- double theta_proton;
- double phi_proton;
-
- double e_fragment;
- double theta_fragment;
- double phi_fragment;
-
- double xtarget = 0;
- double ytarget = 0;
- double CATSX_diff = 0;
- double CATSY_diff = 0;
-
- if (CATS1_X > -1500 && CATS2_X > -1500 && Xf > -1500) {
- CATSX_diff = -(CATS2_X - CATS1_X);
- xtarget = -Xf;
- }
- if (CATS1_Y > -1500 && CATS2_X > -1500 && Yf > -1500) {
- ytarget = Yf;
- CATSY_diff = CATS2_X - CATS1_X;
- }
+ double e_alpha = -1000;
+ double theta_alpha = -1000;
+ double phi_alpha = -1000;
- TVector3 BeamImpact(xtarget, ytarget, 0);
- TVector3 BeamDirection(CATSX_diff, CATSY_diff, 497.1);
+ double e_li = -1000;
+ double theta_li = -1000;
+ double phi_li = -1000;
if (M2->Si_E.size() == 2) {
- if ((cuta->IsInside(M2->CsI_E[0], M2->Si_E[0]) && cutp->IsInside(M2->CsI_E[1], M2->Si_E[1]))) {
+ if ((cuta->IsInside(M2->CsI_E[0], M2->Si_E[0])) && (cutli->IsInside(M2->CsI_E[1], M2->Si_E[1]))) {
// Particle1 = alpha
- TVector3 PositionOfInteraction_alpha(M2_X->at(0), M2_Y->at(0), M2_Z->at(0));
- TVector3 HitDirection_alpha = PositionOfInteraction_alpha - BeamImpact;
e_alpha = ELab->at(0);
- theta_alpha = HitDirection_alpha.Angle(BeamDirection);
- phi_alpha = HitDirection_alpha.Phi();
-
- // Particle2 = proton
- TVector3 PositionOfInteraction_proton(M2_X->at(1), M2_Y->at(1), M2_Z->at(1));
- TVector3 HitDirection_proton = PositionOfInteraction_proton - BeamImpact;
- e_proton = ELab->at(1);
- theta_proton = HitDirection_proton.Angle(BeamDirection);
- phi_proton = HitDirection_proton.Phi();
+ theta_alpha = ThetaLab->at(0) * M_PI / 180.;
+ phi_alpha = PhiLab->at(0);
+ if (phi_alpha + 180. > 360)
+ phi_alpha = phi_alpha - 180;
+ else
+ phi_alpha = phi_alpha + 180;
+ phi_alpha = phi_alpha * M_PI / 180.;
+
+ // Particle2 = 7li
+ e_li = ELab->at(1);
+ theta_li = ThetaLab->at(1) * M_PI / 180.;
+ phi_li = PhiLab->at(1);
+ if (phi_li + 180. > 360)
+ phi_li = phi_li - 180;
+ else
+ phi_li = phi_li + 180;
+ phi_li = phi_li * M_PI / 180.;
+ hExMM->Fill(Ex->at(0));
}
- else if (cuta->IsInside(M2->CsI_E[1], M2->Si_E[1]) && cutp->IsInside(M2->CsI_E[0], M2->Si_E[0])) {
+ else if ((cuta->IsInside(M2->CsI_E[1], M2->Si_E[1])) && (cutli->IsInside(M2->CsI_E[0], M2->Si_E[0]))) {
+ hExMM->Fill(Ex->at(1));
// Particle2 = alpha
- TVector3 PositionOfInteraction_alpha(M2_X->at(1), M2_Y->at(1), M2_Z->at(1));
- TVector3 HitDirection_alpha = PositionOfInteraction_alpha - BeamImpact;
e_alpha = ELab->at(1);
- theta_alpha = HitDirection_alpha.Angle(BeamDirection);
- phi_alpha = HitDirection_alpha.Phi();
-
- // Particle1 = proton
- TVector3 PositionOfInteraction_proton(M2_X->at(0), M2_Y->at(0), M2_Z->at(0));
- TVector3 HitDirection_proton = PositionOfInteraction_proton - BeamImpact;
- e_proton = ELab->at(0);
- theta_proton = HitDirection_proton.Angle(BeamDirection);
- phi_proton = HitDirection_proton.Phi();
+ theta_alpha = ThetaLab->at(1) * M_PI / 180.;
+ phi_alpha = PhiLab->at(1);
+ if (phi_alpha + 180. > 360)
+ phi_alpha = phi_alpha - 180;
+ else
+ phi_alpha = phi_alpha + 180;
+ phi_alpha = phi_alpha * M_PI / 180.;
+
+ // Particle1 = 7li
+ e_li = ELab->at(0);
+ theta_li = ThetaLab->at(0) * M_PI / 180.;
+ phi_li = PhiLab->at(0);
+ if (phi_li + 180. > 360)
+ phi_li = phi_li - 180;
+ else
+ phi_li = phi_li + 180;
+ phi_li = phi_li * M_PI / 180.;
}
- ////////////////////////////////////////////////////////////////////////////////
- double gamma_alpha = e_alpha / m_alpha + 1;
- double beta_alpha = sqrt(1 - 1 / (pow(gamma_alpha, 2.)));
- double p_alpha_mag = gamma_alpha * m_alpha * beta_alpha;
- double etot_alpha = sqrt(pow(p_alpha_mag, 2) + pow(m_alpha, 2));
-
- TVector3 p_alpha(sin(theta_alpha) * cos(phi_alpha), // px
- sin(theta_alpha) * sin(phi_alpha), // py
- cos(theta_alpha)); // pz
- p_alpha.SetMag(p_alpha_mag);
- LV_alpha.SetPxPyPzE(p_alpha.x(), p_alpha.y(), p_alpha.z(), etot_alpha);
-
- ////////////////////////////////////////////////////////////////////////////////
- double gamma_proton = e_proton / m_proton + 1;
- double beta_proton = sqrt(1 - 1 / (pow(gamma_proton, 2.)));
- double p_proton_mag = gamma_proton * m_proton * beta_proton;
- double etot_proton = sqrt(pow(p_proton_mag, 2) + pow(m_proton, 2));
-
- TVector3 p_proton(sin(theta_proton) * cos(phi_proton), // px
- sin(theta_proton) * sin(phi_proton), // py
- cos(theta_proton)); // pz
- p_proton.SetMag(p_proton_mag);
- LV_proton.SetPxPyPzE(p_proton.x(), p_proton.y(), p_proton.z(), etot_proton);
- TLorentzVector LV_total = LV_alpha + LV_proton;
-
- double ExTot = LV_total.Mag() - m_alpha - m_proton;
- hExTot->Fill(ExTot);
- if (Ex.size() >= 2) {
- hEx->Fill(Ex.at(0));
- hEx->Fill(Ex.at(1));
+ if (e_alpha > 0 && e_li > 0) {
+ ////////////////////////////////////////////////////////////////////////////////
+ double gamma_alpha = e_alpha / m_alpha + 1;
+ double beta_alpha = sqrt(1 - 1 / (pow(gamma_alpha, 2.)));
+ cout << i << endl;
+ cout << "alpha: " << e_alpha << " " << beta_alpha << " " << theta_alpha << endl;
+ double p_alpha_mag = gamma_alpha * m_alpha * beta_alpha;
+ double etot_alpha = sqrt(pow(p_alpha_mag, 2) + pow(m_alpha, 2));
+
+ TVector3 p_alpha(sin(theta_alpha) * cos(phi_alpha), // px
+ sin(theta_alpha) * sin(phi_alpha), // py
+ cos(theta_alpha)); // pz
+ p_alpha.SetMag(p_alpha_mag);
+ LV_alpha.SetPxPyPzE(p_alpha.x(), p_alpha.y(), p_alpha.z(), etot_alpha);
+
+ // double etot_alpha = m_alpha + e_alpha;
+ // double p_alpha_mag = sqrt(etot_alpha * etot_alpha - m_alpha * m_alpha);
+ // LV_alpha.SetPxPyPzE(p_alpha_mag * sin(theta_alpha) * cos(phi_alpha), p_alpha_mag * sin(phi_alpha),
+ // p_alpha_mag * cos(theta_alpha), etot_alpha);
+
+ ////////////////////////////////////////////////////////////////////////////////
+ double gamma_li = e_li / m_7Li + 1;
+ double beta_li = sqrt(1 - 1 / (pow(gamma_li, 2.)));
+ cout << "li: " << e_li << " " << beta_li << " " << theta_li << endl;
+ double p_li_mag = gamma_li * m_7Li * beta_li;
+ double etot_li = sqrt(pow(p_li_mag, 2) + pow(m_7Li, 2));
+
+ TVector3 p_li(sin(theta_li) * cos(phi_li), // px
+ sin(theta_li) * sin(phi_li), // py
+ cos(theta_li)); // pz
+ p_li.SetMag(p_li_mag);
+ LV_7Li.SetPxPyPzE(p_li.x(), p_li.y(), p_li.z(), etot_li);
+
+ // double etot_li = m_7Li + e_li;
+ // double p_li_mag = sqrt(etot_li * etot_li - m_7Li * m_7Li);
+ // LV_7Li.SetPxPyPzE(p_li_mag * sin(theta_li) * cos(phi_li), p_li_mag * sin(phi_li), p_li_mag *
+ // cos(theta_li),
+ // etot_li);
+
+ LV_p.SetPxPyPzE(0, 0, 0, mp);
+ TLorentzVector LV_total = LV_alpha + LV_7Li - LV_p;
+
+ double ExTot = LV_total.Mag() - m_10Be;
+ cout << ExTot << endl;
+ hExInv->Fill(ExTot);
+ hE1E2->Fill(e_alpha, e_li);
+ hELabThetaLab->Fill(theta_alpha * 180 / M_PI, e_alpha);
+ hELabThetaLab->Fill(theta_li * 180 / M_PI, e_li);
}
- hE1E2->Fill(e_alpha, e_proton);
}
}
TCanvas* c1 = new TCanvas();
- hExTot->Draw();
- hEx->Draw("same");
+ gStyle->SetOptStat(0);
+ hExInv->Draw();
+ hExMM->SetLineColor(kRed);
+ hExMM->Draw("same");
+
+ TF1* fitExTot = new TF1("fitExTot", "gaus", hExInv->GetXaxis()->GetXmin(), hExInv->GetXaxis()->GetXmax());
+ hExInv->Fit(fitExTot, "Q");
+ double meanExTot = fitExTot->GetParameter(1);
+ std::cout << "Center of hExInv: " << meanExTot << " MeV" << std::endl;
+
+ double fitRangeMin = 0.0;
+ double fitRangeMax = 15.0;
+ TF1* fitEx = new TF1("fitEx", "gaus", fitRangeMin, fitRangeMax);
+ hExMM->Fit(fitEx, "", "", -10, 10);
+ double meanEx = fitEx->GetParameter(1);
+ std::cout << "Center of hEx: " << meanEx << " MeV" << std::endl;
+
+ TLegend* legend = new TLegend(0.75, 0.75, 0.85, 0.85);
+ legend->AddEntry(hExInv, "hExInv", "l"); // "l" for line
+ legend->AddEntry(hExMM, "hExMM", "l");
+ legend->Draw();
+
+ TCanvas* c2 = new TCanvas();
+ hELabThetaLab->Draw("col");
}
-~
diff --git a/Projects/e870/reaction/10Bepalpha.reaction b/Projects/e870/reaction/10Bepalpha.reaction
index a4b5f079592b21fa8e3deaba9be601711f87b4f3..3945abf5819ef29662dfec91a4a517eb9ac3a662 100644
--- a/Projects/e870/reaction/10Bepalpha.reaction
+++ b/Projects/e870/reaction/10Bepalpha.reaction
@@ -13,9 +13,9 @@ Beam
SigmaX= 6.216 mm
SigmaY= 6.109 mm
% Playground:
- %SigmaEnergy= ?? MeV
- %SigmaX= ?? mm
- %SigmaY= ?? mm
+ %SigmaEnergy= 23.2 MeV
+ %SigmaX= 7 mm
+ %SigmaY= 6.8 mm
MeanThetaX= 0 deg
MeanPhiY= 0 deg