From dbd71bd49a89a603e125742d69161d70655b1cea Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?GIRARD=20ALCINDOR=20Val=C3=A9rian?=
<girard-alcindor@ijclab.in2p3.fr>
Date: Thu, 25 Jan 2024 23:12:16 +0100
Subject: [PATCH] Update e870
---
Projects/e870/Analysis.cxx | 4 +
Projects/e870/Analysis.h | 1 +
Projects/e870/macros/DrawExInvariant.cxx | 231 ++++++++++++-----------
3 files changed, 121 insertions(+), 115 deletions(-)
diff --git a/Projects/e870/Analysis.cxx b/Projects/e870/Analysis.cxx
index 928df1b41..a3345f779 100644
--- a/Projects/e870/Analysis.cxx
+++ b/Projects/e870/Analysis.cxx
@@ -134,6 +134,7 @@ void Analysis::TreatEvent() {
TVector3 HitDirection = M2->GetPositionOfInteraction(countMust2) - BeamImpact;
TVector3 Coords = M2->GetPositionOfInteraction(countMust2);
double Theta = HitDirection.Angle(BeamDirection);
+ double Phi = HitDirection.Phi();
X.push_back(Coords.x());
Y.push_back(Coords.y());
Z.push_back(Coords.z());
@@ -170,6 +171,7 @@ void Analysis::TreatEvent() {
// 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);
@@ -186,6 +188,7 @@ 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"); Like before
@@ -237,6 +240,7 @@ void Analysis::ReInitValue() {
dE = -1000;
ELab.clear();
ThetaLab.clear();
+ PhiLab.clear();
Ex.clear();
}
diff --git a/Projects/e870/Analysis.h b/Projects/e870/Analysis.h
index 6ba81025d..afe2d339f 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;
diff --git a/Projects/e870/macros/DrawExInvariant.cxx b/Projects/e870/macros/DrawExInvariant.cxx
index 0d2aa1f11..c624a907a 100644
--- a/Projects/e870/macros/DrawExInvariant.cxx
+++ b/Projects/e870/macros/DrawExInvariant.cxx
@@ -1,16 +1,15 @@
void DrawExInvariant() {
////////////////////////////////////////////////////////////////////////////////
- TLorentzVector LV_alpha, LV_7Li;
- // TLorentzVector LV_alpha, LV_proton, LV_fragment;
- 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_proton = mu + 7.289; // Mev/c2
- // double m_fragment = 44 * mu - 3.755; // 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");
- t->Add("../../Outputs/Analysis/testnew.root"); // Name of simulated file
+ t->Add("../../Outputs/Analysis/Ana10Bepalpha_38AMeV_CH2.root"); // Name of simulated file
// YOUR ALPHA AND Li CUT
TFile* fcuta = new TFile("cuts/CUT_alpha.root", "read");
@@ -22,19 +21,16 @@ void DrawExInvariant() {
TMust2Physics* M2 = new TMust2Physics();
std::vector<double>* ELab = 0;
std::vector<double>* ThetaLab = 0;
- std::vector<double>* X = 0;
- std::vector<double>* Y = 0;
- std::vector<double>* Z = 0;
+ std::vector<double>* PhiLab = 0;
std::vector<double>* Ex = 0;
TBranch* b_ELab;
TBranch* b_ThetaLab;
+ TBranch* b_PhiLab;
TBranch* b_X;
TBranch* b_Y;
TBranch* b_Z;
TBranch* b_Ex;
- // float CATS1_X, CATS1_Y, CATS2_X, CATS2_Y, Xf, Yf;
-
t->SetBranchStatus("*", false);
t->SetBranchStatus("MUST2", true);
t->SetBranchAddress("MUST2", &M2);
@@ -42,145 +38,150 @@ void DrawExInvariant() {
t->SetBranchAddress("ELab", &ELab, &b_ELab);
t->SetBranchStatus("ThetaLab", true);
t->SetBranchAddress("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("X", "true");
- t->SetBranchAddress("X", &X, &b_X);
- t->SetBranchStatus("Y", "true");
- t->SetBranchAddress("Y", &Y, &b_Y);
- t->SetBranchStatus("Z", "true");
- t->SetBranchAddress("Z", &Z, &b_Z);
+ t->SetBranchStatus("PhiLab", true);
+ t->SetBranchAddress("PhiLab", &PhiLab, &b_PhiLab);
t->SetBranchStatus("Ex", "true");
t->SetBranchAddress("Ex", &Ex, &b_Ex);
- TH1F* hExTot = new TH1F("hExTot", "hExTot", 1000, -100, 100);
- TH1F* hEx = new TH1F("hEx", "hEx", 1000, -100, 100);
+ 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_li;
- double theta_li;
- double phi_li;
-
- double xtarget = 0;
- double ytarget = 0;
- // double CATSX_diff = 0;
- // double CATSY_diff = 0;
+ double e_alpha = -1000;
+ double theta_alpha = -1000;
+ double phi_alpha = -1000;
- // 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(0,0,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]) && cutli->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(X->at(0), Y->at(0), 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();
+ 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
- TVector3 PositionOfInteraction_li(X->at(1), Y->at(1), Z->at(1));
- TVector3 HitDirection_li = PositionOfInteraction_li - BeamImpact;
e_li = ELab->at(1);
- theta_li = HitDirection_li.Angle(BeamDirection);
- phi_li = HitDirection_li.Phi();
+ 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]) && cutli->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(X->at(1), Y->at(1), 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();
+ 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
- TVector3 PositionOfInteraction_li(X->at(0), Y->at(0), Z->at(0));
- TVector3 HitDirection_li = PositionOfInteraction_li - BeamImpact;
e_li = ELab->at(0);
- theta_li = HitDirection_li.Angle(BeamDirection);
- phi_li = HitDirection_li.Phi();
+ 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_li = e_li / m_7Li + 1;
- double beta_li = sqrt(1 - 1 / (pow(gamma_li, 2.)));
- 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);
- TLorentzVector LV_total = LV_alpha + LV_7Li;
-
- double ExTot = LV_total.Mag() - m_alpha - m_7Li;
- 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_li);
}
}
TCanvas* c1 = new TCanvas();
gStyle->SetOptStat(0);
- hExTot->Draw();
- hEx->SetLineColor(kRed);
- hEx->Draw("same");
+ hExInv->Draw();
+ hExMM->SetLineColor(kRed);
+ hExMM->Draw("same");
- TF1* fitExTot = new TF1("fitExTot", "gaus", hExTot->GetXaxis()->GetXmin(), hExTot->GetXaxis()->GetXmax());
- hExTot->Fit(fitExTot, "Q");
+ 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 hExTot: " << meanExTot << " MeV" << std::endl;
+ std::cout << "Center of hExInv: " << meanExTot << " MeV" << std::endl;
- double fitRangeMin = 0.0;
- double fitRangeMax = 15.0;
+ double fitRangeMin = 0.0;
+ double fitRangeMax = 15.0;
TF1* fitEx = new TF1("fitEx", "gaus", fitRangeMin, fitRangeMax);
- hEx->Fit(fitEx);
+ 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(hExTot, "hExTot", "l"); // "l" for line
- legend->AddEntry(hEx, "hex", "l");
+ legend->AddEntry(hExInv, "hExInv", "l"); // "l" for line
+ legend->AddEntry(hExMM, "hExMM", "l");
legend->Draw();
-}
+ TCanvas* c2 = new TCanvas();
+ hELabThetaLab->Draw("col");
+}
--
GitLab