diff --git a/NPLib/scripts/build_dict.sh.in b/NPLib/scripts/build_dict.sh.in
index 06950b3ef1abec02aad5f49ffdc32d874bec2c50..e032afeff9cffb70e6746e129e76136542b319e4 100755
--- a/NPLib/scripts/build_dict.sh.in
+++ b/NPLib/scripts/build_dict.sh.in
@@ -38,7 +38,8 @@ if [ $is_root != "" ];
fi
# Get the Major version
-version_major="${version%.*}"
+version_major_string="${version%.*}"
+version_major=$(($version_major_string))
# if before version 4, exit
if [ $version_major -lt 5 ]
diff --git a/NPSimulation/Detectors/MUST2/MUST2Array.cc b/NPSimulation/Detectors/MUST2/MUST2Array.cc
index 831edab0fd9535d954542613cfdf110a5642eb6a..ec47e723258b3c572ef3c83058cbdc7108e8f6a3 100644
--- a/NPSimulation/Detectors/MUST2/MUST2Array.cc
+++ b/NPSimulation/Detectors/MUST2/MUST2Array.cc
@@ -933,7 +933,7 @@ void MUST2Array::ReadSensitive(const G4Event*) {
vector<unsigned int> level = CsIScorer->GetLevel(i);
if (ECsI > ThresholdCsI) {
if (m_CsIOffset[level[0]] == 1)
- m_Event->SetCsIE(level[0], level[1], NPL::EnergyToADC(ECsI, 0, 500, 0, 16384));
+ m_Event->SetCsIE(level[0], level[1], NPL::EnergyToADC(ECsI, 0, 0, 0, 16384));
else
m_Event->SetCsIE(level[0], level[1], NPL::EnergyToADC(ECsI, 0, 250, 8192, 16384));
double timeCsI = RandGauss::shoot(CsIScorer->GetTime(i), ResoTimeMust);
diff --git a/Projects/E805/DetectorConfiguration/MUST2_E805.detector b/Projects/E805/DetectorConfiguration/MUST2_E805.detector
index 98663a5d77eea74774214e348eeb49c2912a101f..0609faabff7d887a6fe1924e5aaa9c28732a8259 100644
--- a/Projects/E805/DetectorConfiguration/MUST2_E805.detector
+++ b/Projects/E805/DetectorConfiguration/MUST2_E805.detector
@@ -1,6 +1,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Target
- THICKNESS= 53.5 micrometer
+ %THICKNESS= 53.5 micrometer
+ THICKNESS= 0.000001 micrometer
ANGLE= 0 deg
RADIUS= 15 mm
MATERIAL= CH2
diff --git a/Projects/E805/project.config b/Projects/E805/project.config
deleted file mode 100644
index f4d49abf59994412aa11998b16b6959b5abd9a3e..0000000000000000000000000000000000000000
--- a/Projects/E805/project.config
+++ /dev/null
@@ -1,8 +0,0 @@
-Project E805
- AnalysisOutput= ./data/NPRoot/Analysis/
- SimulationOutput= ./data/NPRoot/Simulation/
- EnergyLoss= ../../Inputs/EnergyLoss/
- CalibrationOutput= ./data/NPRoot/Calibration/
- Cuts= ./data/NPRoot/Cuts/
-
-
diff --git a/Projects/E805_simu/48Cr_pd_test.reaction b/Projects/E805_simu/48Cr_pd_test.reaction
index ffeab1b9ffc69ddbec7d91f1cb4009eb792023f4..99ae3a287f0f802851376d23328908055a7e3743 100644
--- a/Projects/E805_simu/48Cr_pd_test.reaction
+++ b/Projects/E805_simu/48Cr_pd_test.reaction
@@ -4,12 +4,22 @@
Beam
Particle= 48Cr
Energy= 1511 MeV
- SigmaEnergy= 42 MeV
+ %SigmaEnergy= 42 MeV
+ %ExcitationEnergy= 0 MeV
+ %SigmaThetaX= 0.1 mrad
+ %SigmaPhiY= 0.1 mrad
+ %SigmaX= 5 mm
+ %SigmaY= 5 mm
+ %MeanThetaX= 0 mrad
+ %MeanPhiY= 0 mrad
+ %MeanX= 0 mm
+ %MeanY= 0 mm
+ SigmaEnergy= 0 MeV
ExcitationEnergy= 0 MeV
- SigmaThetaX= 0.1 mrad
- SigmaPhiY= 0.1 mrad
- SigmaX= 5 mm
- SigmaY= 5 mm
+ SigmaThetaX= 0 mrad
+ SigmaPhiY= 0 mrad
+ SigmaX= 0 mm
+ SigmaY= 0 mm
MeanThetaX= 0 mrad
MeanPhiY= 0 mrad
MeanX= 0 mm
@@ -24,7 +34,7 @@ TwoBodyReaction
Target= 1H
Light= 2H
Heavy= 47Cr
- ExcitationEnergyLight= 0.47 MeV
+ ExcitationEnergyLight= 0.0 MeV
ExcitationEnergyHeavy= 0.0 MeV
%CrossSectionPath= CrossSections/48Cr_tot_0.txt CS48Cr
CrossSectionPath= flat.txt CS
diff --git a/Projects/E805_simu/Analysis.cxx b/Projects/E805_simu/Analysis.cxx
index 95867f102122a4e792b2ca886b2706e90c229eb7..cd83a6b1d1de656813ec48f3d266fca08d4c0dbf 100755
--- a/Projects/E805_simu/Analysis.cxx
+++ b/Projects/E805_simu/Analysis.cxx
@@ -43,6 +43,7 @@ void Analysis::Init(){
string beam= NPL::ChangeNameToG4Standard(reaction->GetNucleus1()->GetName());
string heavy_ejectile= NPL::ChangeNameToG4Standard(reaction->GetNucleus4()->GetName());
string light=NPL::ChangeNameToG4Standard(reaction->GetNucleus3()->GetName());
+ std::cout << light << " " << heavy_ejectile << std::endl;
string Reaction_pd_s = "48Cr(p,d)47Cr@1620";
string Reaction_pt_s = "48Cr(p,t)46Cr@1620";
@@ -132,8 +133,11 @@ void Analysis::TreatEvent(){
}
else
Energy = -1000;
- if(Energy > 0)
- Energy = LightTarget["deuteron"].EvaluateInitialEnergy(Energy,TargetThickness*0.5, ThetaNormalTarget);
+ if(Energy > 0){
+ // Energy = LightTarget["deuteron"].EvaluateInitialEnergy(Energy,TargetThickness*0.5, ThetaNormalTarget);
+ Energy = LightTarget["alpha"].EvaluateInitialEnergy(Energy,TargetThickness*0.5, ThetaNormalTarget);
+ }
+
// Evaluate energy using the thickness
diff --git a/Projects/E805_simu/MUST2_E805.detector b/Projects/E805_simu/MUST2_E805.detector
index b15b96f6da402074f8006b7a710ec8a128a79437..1437b43cda0060c66a573e36b31faba4111527d4 100644
--- a/Projects/E805_simu/MUST2_E805.detector
+++ b/Projects/E805_simu/MUST2_E805.detector
@@ -1,6 +1,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Target
- THICKNESS= 53.5 micrometer
+ %THICKNESS= 53.5 micrometer
+ THICKNESS= 0.000001 micrometer
ANGLE= 0 deg
RADIUS= 15 mm
MATERIAL= CH2
diff --git a/Projects/E805_simu/project.config b/Projects/E805_simu/project.config
deleted file mode 100644
index 0edda1e97a786835c74354cfe800128eb48e8261..0000000000000000000000000000000000000000
--- a/Projects/E805_simu/project.config
+++ /dev/null
@@ -1,8 +0,0 @@
-Project E805_simu
- AnalysisOutput= ../E805/data/NPRoot/Analysis/
- SimulationOutput= ../E805/data/NPRoot/Simulation/
- EnergyLoss= ../../Inputs/EnergyLoss/
- CalibrationOutput= ../E805/data/NPRoot/Calibration/
- Cuts= ../E805/data/NPRoot/Cuts/
-
-
diff --git a/Projects/e870/.DS_Store b/Projects/e870/.DS_Store
index f5ddf34afd94504d47c0de268f5e147d9a740da6..c4bac702c96a1b1e9866a310b007dec4c64658b1 100644
Binary files a/Projects/e870/.DS_Store and b/Projects/e870/.DS_Store differ
diff --git a/Projects/e870/Analysis.cxx b/Projects/e870/Analysis.cxx
index db027a980e60e9255dd3a01b51a7321385f2d95f..af55fe7ba7ff505fa1267fe2f1d7178ab39cf748 100644
--- a/Projects/e870/Analysis.cxx
+++ b/Projects/e870/Analysis.cxx
@@ -50,17 +50,18 @@ void Analysis::Init() {
Initial = new TInitialConditions();
ReactionConditions = new TReactionConditions();
}
-
+ // ExogamData = new TExogamData();
InitOutputBranch();
InitInputBranch();
// get MUST2 and Gaspard objects
M2 = (TMust2Physics*)m_DetectorManager->GetDetector("M2Telescope");
- if (!simulation)
- CATS = (TCATSPhysics*)m_DetectorManager->GetDetector("CATSDetector");
+ // if (!simulation)
+ // CATS = (TCATSPhysics*)m_DetectorManager->GetDetector("CATSDetector");
// get reaction information
reaction.ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
OriginalBeamEnergy = reaction.GetBeamEnergy();
+ std::cout << OriginalBeamEnergy << std::endl;
// target thickness
TargetThickness = m_DetectorManager->GetTargetThickness();
string TargetMaterial = m_DetectorManager->GetTargetMaterial();
@@ -70,6 +71,7 @@ void Analysis::Init() {
string beam = NPL::ChangeNameToG4Standard(reaction.GetNucleus1()->GetName());
cout << light << " " << beam << " " << TargetMaterial << " " << TargetThickness << endl;
LightTarget = NPL::EnergyLoss(light + "_" + TargetMaterial + ".G4table", "G4Table", 100);
+ LightMyl = NPL::EnergyLoss(light + "_Mylar.G4table", "G4Table", 100);
LightAl = NPL::EnergyLoss(light + "_Al.G4table", "G4Table", 100);
LightSi = NPL::EnergyLoss(light + "_Si.G4table", "G4Table", 100);
BeamTarget = NPL::EnergyLoss(beam + "_" + TargetMaterial + ".G4table", "G4Table", 100);
@@ -80,17 +82,6 @@ void Analysis::Init() {
<< "MeV " << endl;
reaction.SetBeamEnergy(FinalBeamEnergy);
- if (WindowsThickness) {
- cout << "Cryogenic target with windows" << endl;
- // BeamWindow = new NPL::EnergyLoss(beam + "_" + WindowsMaterial + ".G4table", "G4Table", 100);
- // LightWindow = new NPL::EnergyLoss(light + "_" + WindowsMaterial + ".G4table", "G4Table", 100);
- }
-
- else {
- BeamWindow = NULL;
- LightWindow = NULL;
- }
-
// initialize various parameters
Rand = TRandom3();
DetectorNumber = 0;
@@ -114,19 +105,12 @@ void Analysis::TreatEvent() {
ReInitValue();
double XTarget, YTarget;
TVector3 BeamDirection;
- if (!simulation) {
- XTarget = CATS->GetPositionOnTarget().X();
- YTarget = CATS->GetPositionOnTarget().Y();
- BeamDirection = CATS->GetBeamDirection();
- }
- else {
- XTarget = 0;
- YTarget = 0;
- BeamDirection = TVector3(0, 0, 1);
- // OriginalELab = ReactionConditions->GetKineticEnergy(0);
- // OriginalThetaLab = ReactionConditions->GetTheta(0);
- // BeamEnergy = ReactionConditions->GetBeamEnergy();
- }
+
+ // For simulation
+ XTarget = 0;
+ YTarget = 0;
+ BeamDirection = TVector3(0, 0, 1);
+
BeamImpact = TVector3(XTarget, YTarget, 0);
// determine beam energy for a randomized interaction point in target
// 1% FWHM randominastion (E/100)/2.35
@@ -136,7 +120,8 @@ void Analysis::TreatEvent() {
////////////////////////////////////////////////////////////////////////////
//////////////////////////////// LOOP on MUST2 ////////////////////////////
////////////////////////////////////////////////////////////////////////////
- for (unsigned int countMust2 = 0; countMust2 < M2->Si_E.size(); countMust2++) {
+ int M2_size = M2->Si_E.size();
+ for (unsigned int countMust2 = 0; countMust2 < M2_size; countMust2++) {
/************************************************/
// Part 0 : Get the usefull Data
// MUST2
@@ -168,9 +153,8 @@ void Analysis::TreatEvent() {
if (CsI_E_M2 > 0) {
// The energy in CsI is calculate form dE/dx Table because
Energy = CsI_E_M2;
- if (simulation) {
- Energy = LightAl.EvaluateInitialEnergy(Energy, 0.4 * micrometer, ThetaM2Surface);
- }
+ Energy = LightMyl.EvaluateInitialEnergy(Energy, 3 * micrometer, ThetaM2Surface);
+ Energy = LightAl.EvaluateInitialEnergy(Energy, 0.4 * micrometer, ThetaM2Surface);
Energy += Si_E_M2;
}
else {
@@ -180,7 +164,7 @@ void Analysis::TreatEvent() {
Energy = LightAl.EvaluateInitialEnergy(Energy, 0.4 * micrometer, ThetaM2Surface);
// Evaluate energy using the thickness
// Target Correction
- Energy = LightTarget.EvaluateInitialEnergy(Energy, TargetThickness * 0.5, ThetaNormalTarget);
+ Energy = LightTarget.EvaluateInitialEnergy(Energy, TargetThickness * 0.5, Theta);
// What is written in the tree
ThetaLab.push_back(Theta / deg);
@@ -195,6 +179,7 @@ void Analysis::TreatEvent() {
void Analysis::End() {}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
+ // RootOutput::getInstance()->GetTree()->Branch("Exogam", &ExogamData);
RootOutput::getInstance()->GetTree()->Branch("Ex", &Ex);
RootOutput::getInstance()->GetTree()->Branch("ELab", &ELab);
RootOutput::getInstance()->GetTree()->Branch("ThetaLab", &ThetaLab);
@@ -226,6 +211,10 @@ void Analysis::InitInputBranch() {
RootInput::getInstance()->GetChain()->SetBranchStatus("fRC_*", true);
RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions", &ReactionConditions);
}
+ // RootInput::getInstance()->GetChain()->SetBranchStatus("Exogam", true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("fExo*", true);
+ // RootInput::getInstance()->GetChain()->SetBranchAddress("Exogam", &ExogamData);
+ // RootInput::getInstance()->GetChain()->SetBranchStatus("ZDD", true);
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue() {
diff --git a/Projects/e870/Analysis.h b/Projects/e870/Analysis.h
index 450e560cbbfbdb7f318a4c2194bad69af6204108..753d9fcdbfb81f25e6efff8ba3727f67561884b3 100644
--- a/Projects/e870/Analysis.h
+++ b/Projects/e870/Analysis.h
@@ -30,7 +30,8 @@
#include "TReactionConditions.h"
#include "TMust2Physics.h"
#include "TMugastPhysics.h"
-#include "TCATSPhysics.h"
+// #include "TCATSPhysics.h"
+#include "TExogamData.h"
#include <TRandom3.h>
#include <TVector3.h>
#include <TMath.h>
@@ -64,6 +65,7 @@ class Analysis: public NPL::VAnalysis{
NPL::Reaction reaction;
// Energy loss table: the G4Table are generated by the simulation
NPL::EnergyLoss LightTarget;
+ NPL::EnergyLoss LightMyl;
NPL::EnergyLoss LightAl;
NPL::EnergyLoss LightSi;
NPL::EnergyLoss BeamTarget;
@@ -118,9 +120,11 @@ class Analysis: public NPL::VAnalysis{
// Branches and detectors
TMust2Physics* M2;
TMugastPhysics* MG;
- TCATSPhysics* CATS;
+ // TCATSPhysics* CATS;
bool simulation;
TInitialConditions* Initial;
TReactionConditions* ReactionConditions;
+ TExogamData* ExogamData;
+
};
#endif
diff --git a/Projects/e870/DetectorConfiguration/MUGAST_LISE.detector b/Projects/e870/DetectorConfiguration/MUGAST_LISE.detector
index 670c6f6a727af047791c4db3f01708df155e76b7..a47e82527e81b493cf82ff6f858cd6e155a9c60e 100644
--- a/Projects/e870/DetectorConfiguration/MUGAST_LISE.detector
+++ b/Projects/e870/DetectorConfiguration/MUGAST_LISE.detector
@@ -2,6 +2,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
Target
THICKNESS= 26 micrometer
+ %THICKNESS= 0.00001 micrometer
ANGLE= 0 deg
RADIUS= 13 mm
MATERIAL= CH2
diff --git a/Projects/e870/DetectorConfiguration/MUGAST_LISE_CD2.detector b/Projects/e870/DetectorConfiguration/MUGAST_LISE_CD2.detector
index 020abacbd17a7139fd87ac009d6e6fa893c11ee5..c8cb44e76f2c0332331696758d46f3e691f524ee 100644
--- a/Projects/e870/DetectorConfiguration/MUGAST_LISE_CD2.detector
+++ b/Projects/e870/DetectorConfiguration/MUGAST_LISE_CD2.detector
@@ -36,6 +36,8 @@ M2Telescope
SILI= 0
CSI= 1
VIS= all
+ CsIOffset= 1
+
%%%%%%% Telescope 2 %%%%%%%
M2Telescope
X1_Y1= -114.9 9.68 291.84 mm
@@ -46,6 +48,7 @@ M2Telescope
SILI= 0
CSI= 1
VIS= all
+ CsIOffset= 1
%%%%%%% Telescope 3 %%%%%%%
M2Telescope
@@ -57,18 +60,21 @@ M2Telescope
SILI= 0
CSI= 1
VIS= all
+ CsIOffset= 1
+
+%%%%%%% Telescope 4 %%%%%%%
+M2Telescope
+ X1_Y1= 113.56 -13.18 292.11 mm
+ X1_Y128= 23.23 -13.37 328.15 mm
+ X128_Y1= 102.39 -105.49 263.59 mm
+ X128_Y128= 12.04 -105.69 299.63 mm
+ SI= 1
+ SILI= 0
+ CSI= 1
+ VIS= all
+ CsIOffset= 1
-% %%%%%%% Telescope 4 %%%%%%%
-% M2Telescope
-% X1_Y1= 113.56 -13.18 292.11 mm
-% X1_Y128= 23.23 -13.37 328.15 mm
-% X128_Y1= 102.39 -105.49 263.59 mm
-% X128_Y128= 12.04 -105.69 299.63 mm
-% SI= 1
-% SILI= 0
-% CSI= 1
-% VIS= all
%
% %%%%%%% Telescope 1 %%%%%%%
% M2Telescope
diff --git a/Projects/e870/launch_simu.sh b/Projects/e870/launch_simu.sh
index 9005ba4a2f61168db03b0c9e79258597071a5b67..da74d29302d83ba6a884c2c0d32d0d501b8cd07d 100755
--- a/Projects/e870/launch_simu.sh
+++ b/Projects/e870/launch_simu.sh
@@ -47,9 +47,9 @@
# All of the above:
npcompilation
-npsimulation -D DetectorConfiguration/MUGAST_LISE.detector -E reaction/10Bepalpha.reaction -O test -B run.mac
+npsimulation -D DetectorConfiguration/MUGAST_LISE.detector -E reaction/10Bed6Li.reaction -O test -B run.mac
-npanalysis -D DetectorConfiguration/MUGAST_LISE.detector -E reaction/10Bepalpha.reaction -T $NPTOOL/Outputs/Simulation/test.root SimulatedTree -O test
+npanalysis -D DetectorConfiguration/MUGAST_LISE.detector -E reaction/10Bed6Li.reaction -T $NPTOOL/Outputs/Simulation/test.root SimulatedTree -O test
diff --git a/Projects/e870/macros/DrawExInvariant.cxx b/Projects/e870/macros/DrawExInvariant.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0d191812285a13367b9dcfdda01464145e69b49d
--- /dev/null
+++ b/Projects/e870/macros/DrawExInvariant.cxx
@@ -0,0 +1,166 @@
+void DrawDalitz() {
+ ////////////////////////////////////////////////////////////////////////////////
+ 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
+
+ TChain* t = new TChain("PhysicsTree"); // SimulatedTree
+ t->Add("NPRootA/NPr0371_*a.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");
+
+ ////////////////////////////////////////////////////////////////////////////////
+ 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;
+ 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;
+
+ 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);
+
+ TH2F* hE1E2 = new TH2F("hE1E2", "hE1E2", 1000, 0, 500, 1000, 0, 500);
+
+ ////////////////////////////////////////////////////////////////////////////////
+ // 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;
+ }
+
+ TVector3 BeamImpact(xtarget, ytarget, 0);
+ TVector3 BeamDirection(CATSX_diff, CATSY_diff, 497.1);
+
+ 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]))) {
+ // 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();
+ }
+ else if (cuta->IsInside(M2->CsI_E[1], M2->Si_E[1]) && cutp->IsInside(M2->CsI_E[0], M2->Si_E[0])) {
+ // 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();
+ }
+
+ ////////////////////////////////////////////////////////////////////////////////
+ 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));
+ }
+ hE1E2->Fill(e_alpha, e_proton);
+ }
+ }
+ TCanvas* c1 = new TCanvas();
+ hExTot->Draw();
+ hEx->Draw("same");
+}
+~
diff --git a/Projects/e870/reaction/10Bed6Li.reaction b/Projects/e870/reaction/10Bed6Li.reaction
index ecf3e3413291709c44c3db79303f3791111ce757..dcf2f7459d8d02300451c34a0ae6c792e9e8bf01 100644
--- a/Projects/e870/reaction/10Bed6Li.reaction
+++ b/Projects/e870/reaction/10Bed6Li.reaction
@@ -1,7 +1,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Beam
Particle= 10Be
- Energy= 300 MeV
+ Energy= 150 MeV
SigmaEnergy= 21.6 MeV
MeanThetaX= 0 deg
MeanPhiY= 0 deg
diff --git a/Projects/e870/reaction/10Bepalpha.reaction b/Projects/e870/reaction/10Bepalpha.reaction
index 5c5bf0441c3f73067dd0c3ff1777594c8c3bcf5d..a4b5f079592b21fa8e3deaba9be601711f87b4f3 100644
--- a/Projects/e870/reaction/10Bepalpha.reaction
+++ b/Projects/e870/reaction/10Bepalpha.reaction
@@ -5,13 +5,13 @@ Beam
% Stuff to try and test:
% Perfect case:
- SigmaEnergy= 0 MeV
- SigmaX= 0 mm
- SigmaY= 0 mm
+ % SigmaEnergy= 0 MeV
+ % SigmaX= 0 mm
+ % SigmaY= 0 mm
% Realistic case?
- %SigmaEnergy= 22.8 MeV
- %SigmaX= 6.216 mm
- %SigmaY= 6.109 mm
+ SigmaEnergy= 22.8 MeV
+ SigmaX= 6.216 mm
+ SigmaY= 6.109 mm
% Playground:
%SigmaEnergy= ?? MeV
%SigmaX= ?? mm
@@ -21,13 +21,13 @@ Beam
MeanPhiY= 0 deg
MeanX= 0 mm
MeanY= 0 mm
- %SigmaThetaX= 0.6138 deg
- %SigmaPhiY= 0.3812 deg
+ SigmaThetaX= 0.6138 deg
+ SigmaPhiY= 0.3812 deg
%SigmaEnergy= 0 MeV
%ExcitationEnergy= 0 MeV
- SigmaThetaX= 0 mrad
- SigmaPhiY= 0 mrad
+ %SigmaThetaX= 0 mrad
+ %SigmaPhiY= 0 mrad
%MeanThetaX= 0 mrad
%MeanPhiY= 0 mrad
%MeanX= 0 mm