Skip to content

Commit

Permalink
Merge pull request #8 from vetlewi/master
Browse files Browse the repository at this point in the history 
Fixed some bugs
  • Loading branch information
@vetlewi
vetlewi authored Jul 10, 2016
2 parents afb5ebb + 9f63460 commit 8640d87
Show file tree
Hide file tree
Showing 9 changed files with 108 additions and 154 deletions.
2 changes: 1 addition & 1 deletion app/app.pro
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,4 +40,4 @@ linux {
QMAKE_CLEAN += $$BUILDDIR/app/*.o \
$$BUILDDIR/app/moc_* \
$$BUILDDIR/app/qrc_* \
$$BUILDDIR/app/ui_*
$$BUILDDIR/app/ui_*
27 changes: 10 additions & 17 deletions src/gui/forms/mainwindow.ui
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@
<rect>
<x>0</x>
<y>0</y>
<width>1489</width>
<height>998</height>
<width>1399</width>
<height>849</height>
</rect>
</property>
<property name="windowTitle">
Expand Down Expand Up @@ -38,6 +38,9 @@
<layout class="QVBoxLayout" name="verticalLayout_4">
<item>
<layout class="QVBoxLayout" name="verticalLayout_3">
<property name="sizeConstraint">
<enum>QLayout::SetFixedSize</enum>
</property>
<item>
<layout class="QHBoxLayout" name="Layout_beamInfo">
<item>
Expand Down Expand Up @@ -553,27 +556,17 @@
<property name="orientation">
<enum>Qt::Vertical</enum>
</property>
<property name="sizeType">
<enum>QSizePolicy::MinimumExpanding</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>198</width>
<height>13</height>
<height>0</height>
</size>
</property>
</spacer>
</item>
<item>
<widget class="QLabel" name="OCLlogo">
<property name="text">
<string/>
</property>
<property name="pixmap">
<pixmap resource="../../../resources/resorces.qrc">:/media/clab-logo-med.gif</pixmap>
</property>
<property name="scaledContents">
<bool>true</bool>
</property>
</widget>
</item>
</layout>
</widget>
</item>
Expand Down Expand Up @@ -620,7 +613,7 @@
<rect>
<x>0</x>
<y>0</y>
<width>1489</width>
<width>1399</width>
<height>22</height>
</rect>
</property>
Expand Down
127 changes: 1 addition & 126 deletions src/gui/src/mainwindow.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -329,136 +329,11 @@ void MainWindow::run()
table.Reset();

double angle = (ui->StripNumbr->value()*2 + 40)*PI/180;
if (ui->BwdAngRButton){
if (ui->BwdAngRButton->isChecked()){
angle = PI - angle;
}

emit operate(angle, ui->protons->isChecked(), ui->deutrons->isChecked(), ui->tritons->isChecked(), ui->He3s->isChecked(), ui->alphas->isChecked());

/*
if (ui->protons->isChecked()){
if (worker->Run(Ex, dE, E, coeff, angle, 1, 1)){
QVector<double> x, y, ex;
for (int i = 0 ; i < dE.size() ; ++i){
if (E[i] >= 0.05){
ex.push_back(Ex[i]);
x.push_back(E[i]);
y.push_back(dE[i]);
}
}
table.setCoeff(coeff, TableMakerHTML::Proton);
QCPCurve *nCurve = makeCurve(x, y, QPen(QColor(255,0,0)), QString("%1%2)").arg(legend).arg("p"));
ui->plotTab->addPlottable(nCurve);
}
if (worker->Known(Ex, dE, E, d_dE, d_E, angle, 1, 1)){
table.setData(Ex, dE, d_dE, E, d_E, TableMakerHTML::Proton);
QCPGraph *nGraph = makeGraph(E, d_E, dE, d_dE, QPen(QColor(255,0,0)));
ui->plotTab->addPlottable(nGraph);
ui->plotTab->plottable()->removeFromLegend();
}
ui->plotTab->replot();
ui->plotTab->show();
}
if (ui->deutrons->isChecked()){
if (worker->Run(Ex, dE, E, coeff, angle, 2, 1)){
QVector<double> x, y, ex;
for (int i = 0 ; i < dE.size() ; ++i){
if (E[i] >= 0.05){
ex.push_back(Ex[i]);
x.push_back(E[i]);
y.push_back(dE[i]);
}
}
table.setCoeff(coeff, TableMakerHTML::Deutron);
QCPCurve *nCurve = makeCurve(x, y, QPen(QColor(0,255,0)), QString("%1%2)").arg(legend).arg("d"));
ui->plotTab->addPlottable(nCurve);
}
if (worker->Known(Ex, dE, E, d_dE, d_E, angle, 2, 1)){
table.setData(Ex, dE, d_dE, E, d_E, TableMakerHTML::Deutron);
QCPGraph *nGraph = makeGraph(E, d_E, dE, d_dE, QPen(QColor(0,255,0)));
ui->plotTab->addPlottable(nGraph);
ui->plotTab->plottable()->removeFromLegend();
}
ui->plotTab->replot();
ui->plotTab->show();
}
if (ui->tritons->isChecked()){
if (worker->Run(Ex, dE, E, coeff, angle, 3, 1)){
QVector<double> x, y, ex;
for (int i = 0 ; i < dE.size() ; ++i){
if (E[i] >= 0.05){
ex.push_back(Ex[i]);
x.push_back(E[i]);
y.push_back(dE[i]);
}
}
table.setCoeff(coeff, TableMakerHTML::Triton);
QCPCurve *nCurve = makeCurve(x, y, QPen(QColor(128,128,128)), QString("%1%2)").arg(legend).arg("t"));
ui->plotTab->addPlottable(nCurve);
}
if (worker->Known(Ex, dE, E, d_dE, d_E, angle, 3, 1)){
table.setData(Ex, dE, d_dE, E, d_E, TableMakerHTML::Triton);
QCPGraph *nGraph = makeGraph(E, d_E, dE, d_dE, QPen(QColor(128,128,128)));
ui->plotTab->addPlottable(nGraph);
ui->plotTab->plottable()->removeFromLegend();
}
ui->plotTab->replot();
ui->plotTab->show();
}
if (ui->He3s->isChecked()){
if (worker->Run(Ex, dE, E, coeff, angle, 3, 2)){
QVector<double> x, y, ex;
for (int i = 0 ; i < dE.size() ; ++i){
if (E[i] >= 0.05){
ex.push_back(Ex[i]);
x.push_back(E[i]);
y.push_back(dE[i]);
}
}
table.setCoeff(coeff, TableMakerHTML::Helium3);
QCPCurve *nCurve = makeCurve(x, y, QPen(QColor(0,0,255)), QString("%1%2)").arg(legend).arg("He3"));
ui->plotTab->addPlottable(nCurve);
}
if (worker->Known(Ex, dE, E, d_dE, d_E, angle, 3, 2)){
table.setData(Ex, dE, d_dE, E, d_E, TableMakerHTML::Helium3);
QCPGraph *nGraph = makeGraph(E, d_E, dE, d_dE, QPen(QColor(0,0,255)));
ui->plotTab->addPlottable(nGraph);
ui->plotTab->plottable()->removeFromLegend();
}
ui->plotTab->replot();
ui->plotTab->show();
}
if (ui->alphas->isChecked()){
if (worker->Run(Ex, dE, E, coeff, angle, 4, 2)){
QVector<double> x, y, ex;
for (int i = 0 ; i < dE.size() ; ++i){
if (E[i] >= 0.05){
ex.push_back(Ex[i]);
x.push_back(E[i]);
y.push_back(dE[i]);
}
}
table.setCoeff(coeff, TableMakerHTML::Alpha);
QCPCurve *nCurve = makeCurve(x, y, QPen(QColor(0,0,0)), QString("%1%2)").arg(legend).arg("α"));
ui->plotTab->addPlottable(nCurve);
}
if (worker->Known(Ex, dE, E, d_dE, d_E, angle, 4, 2)){
table.setData(Ex, dE, d_dE, E, d_E, TableMakerHTML::Alpha);
QCPGraph *nGraph = makeGraph(E, d_E, dE, d_dE, QPen(QColor(0,0,0)));
ui->plotTab->addPlottable(nGraph);
ui->plotTab->plottable()->removeFromLegend();
}
ui->plotTab->replot();
ui->plotTab->show();
}
ui->webView->setHtml(table.getHTMLCode());
ui->plotTab->rescaleAxes();
emit runDialog->Finished();
*/
}

void MainWindow::WorkFinished()
Expand Down
34 changes: 34 additions & 0 deletions src/kinematics/include/RelScatter.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
#ifndef RELSCATTER_H
#define RELSCATTER_H


class Particle;
#include "Scattering.h"

class RelScatter : public Scattering
{
public:
//! Constuctor.
RelScatter(Particle *pA, /*!< Incident particle. */
Particle *pX, /*!< Target particle. */
Particle *pY, /*!< Residual particle. */
Particle *pB /*!< Fragement particle. */);

//! Destructor.
~RelScatter();

//! Evaluate energy of residual particle.
/*! \return Energy of particle Y after scattering.
*/
double EvaluateY(const double &E, /*!< Incident energy. */
const double &theta, /*!< Scattering angle of particle Y */
const double &Ex /*!< Excitation of particle B. */) const;


//! Calculate the maximum excitation energy of particle B at given scattering angle and residual energy.
double FindMaxEx(const double &E, /*!< Incident energy. */
const double &theta /*!< Scattering angle. */) const;

};

#endif // RELSCATTER_H
4 changes: 2 additions & 2 deletions src/kinematics/include/Scattering.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,8 +37,8 @@ class Scattering
//! Function for calculating maximum excitation energy of particle B.
/*! \return The highest posible positive excitation energy for particle B.
*/
//virtual double FindMaxEx(const double &E, /*!< Energy of incident particle A. */
// const double &theta /*!< Outgoing angle of the fragment. */) const=0;
virtual double FindMaxEx(const double &E, /*!< Energy of incident particle A. */
const double &theta /*!< Outgoing angle of the fragment. */) const=0;

// Functions to set particles.
//! Set particle A.
Expand Down
49 changes: 49 additions & 0 deletions src/kinematics/src/RelScatter.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,49 @@
#include "RelScatter.h"

#include "Particle.h"

#include <cmath>


RelScatter::RelScatter(Particle *pA, Particle *pX, Particle *pY, Particle *pB)
: Scattering(pA, pX, pY, pB)
{

}

RelScatter::~RelScatter(){ }

double RelScatter::EvaluateY(const double &E, const double &theta, const double &Ex) const
{
double m1 = X->GetM_MeV(), m2 = A->GetM_MeV(), m3 = Y->GetM_MeV(), m4 = B->GetM_MeV();

double m4ex = m4 + Ex;

double s = (m1+m2)*(m1+m2) + 2*m1*E;

double pcm = sqrt(((s - m1*m1 - m2*m2)*(s - m1*m1 - m2*m2) - 4*m1*m1*m2*m2)/(4*s));

double chi = log((pcm + sqrt(pcm*pcm + m1*m1))/m1);

double pcm2 = sqrt(((s - m3*m3 - m4ex*m4ex)*(s - m3*m3 - m4ex*m4ex) - 4*m3*m3*m4ex*m4ex)/(4*s));

double p3 = sqrt(pcm2*pcm2 + m3*m3)*cos(theta)*sinh(chi)+sqrt(pcm2*pcm2 - pow(m3*sin(theta)*sinh(chi),2))*cosh(chi);
p3 /= 1 + sin(theta)*sin(theta)*sinh(chi)*sinh(chi);

double E3 = sqrt(p3*p3 + m3*m3);
return E3 - m3;
}

double RelScatter::FindMaxEx(const double &E, const double &theta) const
{
double m1 = X->GetM_MeV(), m2 = A->GetM_MeV(), m3 = Y->GetM_MeV(), m4 = B->GetM_MeV();

double s = (m1+m2)*(m1+m2) + 2*m1*E;

double pcm = sqrt(((s - m1*m1 - m2*m2)*(s - m1*m1 - m2*m2) - 4*m1*m1*m2*m2)/(4*s));

double chi = log((pcm + sqrt(pcm*pcm + m1*m1))/m1);

double Ex = sqrt(s - 2*sqrt(s)*m3*sqrt(1 + pow(sin(theta)*sinh(chi), 2)) + m3*m3) - m4;
return Ex;
}
2 changes: 2 additions & 0 deletions src/src.pro
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ SOURCES += gui/src/mainwindow.cpp \
kinematics/src/StoppingPower.cpp \
kinematics/src/Ziegler1985.cpp \
kinematics/src/ZieglerComp.cpp \
kinematics/src/RelScatter.cpp \
math/src/AbstractFunction.cpp \
math/src/Matrix.cpp \
math/src/PolyD2.cpp \
Expand Down Expand Up @@ -79,6 +80,7 @@ HEADERS += gui/include/mainwindow.h \
kinematics/include/StoppingPower.h \
kinematics/include/Ziegler1985.h \
kinematics/include/ZieglerComp.h \
kinematics/include/RelScatter.h \
math/include/AbstractFunction.h \
math/include/Matrix.h \
math/include/PolyD2.h \
Expand Down
7 changes: 4 additions & 3 deletions src/support/src/runsystem.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@

#include "Scattering.h"
#include "DickNorbury.h"
#include "RelScatter.h"
#include "StoppingPower.h"
#include "Ziegler1985.h"
#include "ame2012_masses.h"
Expand Down Expand Up @@ -48,7 +49,7 @@ RunSystem::~RunSystem()

int RunSystem::Run(const double &Energy, const double &Angle) const
{
DickNorbury *scat = new DickNorbury(beam, scatIso, fragment, residual);
RelScatter *scat = new RelScatter(beam, scatIso, fragment, residual);
Ziegler1985 *stopTargetB = new Ziegler1985(target, beam);
Ziegler1985 *stopTargetF = new Ziegler1985(target, fragment);
Ziegler1985 *stopAbsor = new Ziegler1985(absorber, fragment);
Expand All @@ -57,9 +58,9 @@ int RunSystem::Run(const double &Energy, const double &Angle) const

double Ehalf = stopTargetB->Loss(Energy, target->GetWidth(), INTPOINTS);
double Ehole = stopTargetB->Loss(Energy, INTPOINTS);
double Exmax = get_Q_keV(beam->GetA(), beam->GetZ(), scatIso->GetA(), scatIso->GetZ(), fragment->GetA(), fragment->GetZ())/1000;
double Exmax = scat->FindMaxEx(Ehalf, Angle);

Exmax += Energy;
//Exmax += Energy;
double dEX = Exmax/double(POINTS - 1);
double f, m, b;
double df, dm, db;
Expand Down
10 changes: 5 additions & 5 deletions src/support/src/worker.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#include "Scattering.h"
#include "DickNorbury.h"
#include "Iterative.h"
#include "LNScattering.h"
#include "RelScatter.h"

#include "StoppingPower.h"
#include "Ziegler1985.h"
Expand Down Expand Up @@ -112,7 +112,7 @@ bool Worker::Curve(QVector<double> &Ex, QVector<double> &dE, QVector<double> &E,
Material *dEdet = new Material(theTelescope->dEdetector.Z, Get_mm2(theTelescope->dEdetector.Z), theTelescope->dEdetector.width/cos(incAngle), Unit2MatUnit(theTelescope->dEdetector.unit));
Material *Edet = new Material(theTelescope->Edetector.Z, Get_mm2(theTelescope->Edetector.Z), theTelescope->Edetector.width/cos(incAngle), Unit2MatUnit(theTelescope->Edetector.unit));

Scattering *scat = new LNScattering(beam, scatIso, fragment, residual);//new Iterative(beam, scatIso, fragment, residual);
RelScatter *scat = new RelScatter(beam, scatIso, fragment, residual);//new Iterative(beam, scatIso, fragment, residual);

// Setting up stopping power for the target.
StoppingPower *stopTargetB;
Expand Down Expand Up @@ -181,7 +181,7 @@ bool Worker::Curve(QVector<double> &Ex, QVector<double> &dE, QVector<double> &E,
double Ehalf = stopTargetB->Loss(E_beam, target->GetWidth(tUnit)/2., INTPOINTS);
double Ewhole = stopTargetB->Loss(E_beam, INTPOINTS);

double dEx = (Ehalf + get_Q_keV(theBeam->A, theBeam->Z, theTarget->A, theTarget->Z, fA, fZ)/1000.)/double(POINTS - 1);
double dEx = scat->FindMaxEx(Ewhole, Angle)/double(POINTS - 1);

if ((Ehalf + get_Q_keV(theBeam->A, theBeam->Z, theTarget->A, theTarget->Z, fA, fZ)/1000.)<0){
// Particles.
Expand Down Expand Up @@ -364,7 +364,7 @@ bool Worker::Known(QVector<double> &Ex, QVector<double> &dE, QVector<double> &E,
Material *dEdet = new Material(theTelescope->dEdetector.Z, Get_mm2(theTelescope->dEdetector.Z), theTelescope->dEdetector.width/cos(incAngle), Unit2MatUnit(theTelescope->dEdetector.unit));
Material *Edet = new Material(theTelescope->Edetector.Z, Get_mm2(theTelescope->Edetector.Z), theTelescope->Edetector.width/cos(incAngle), Unit2MatUnit(theTelescope->Edetector.unit));

Scattering *scat = new LNScattering(beam, scatIso, fragment, residual);//new Iterative(beam, scatIso, fragment, residual);
RelScatter *scat = new RelScatter(beam, scatIso, fragment, residual);//new Iterative(beam, scatIso, fragment, residual);

// Setting up stopping power for the target.
StoppingPower *stopTargetB;
Expand Down Expand Up @@ -434,7 +434,7 @@ bool Worker::Known(QVector<double> &Ex, QVector<double> &dE, QVector<double> &E,
double Ehalf = stopTargetB->Loss(E_beam, target->GetWidth(tUnit)/2., INTPOINTS);
double Ewhole = stopTargetB->Loss(E_beam, INTPOINTS);

double Exmax = Ehalf + get_Q_keV(theBeam->A, theBeam->Z, theTarget->A, theTarget->Z, fA, fZ)/1000.;
double Exmax = scat->FindMaxEx(Ewhole, Angle);

if ((Ehalf + get_Q_keV(theBeam->A, theBeam->Z, theTarget->A, theTarget->Z, fA, fZ)/1000.)<0){
// Particles.
Expand Down

0 comments on commit 8640d87

Please sign in to comment.