Commit a32ab3d5 authored by ext-calvo_p's avatar ext-calvo_p
Browse files

Resolve "Loss files overwritten for collimators"

parent 266394fa
......@@ -78,11 +78,11 @@ BeamStrippingPhysics::BeamStrippingPhysics(const std::string& name, ElementBase*
bunchToMatStat_m(0),
stoppedPartStat_m(0),
rediffusedStat_m(0),
locPartsInMat_m(0)
totalPartsInMat_m(0)
{
bstp_m = dynamic_cast<BeamStripping* >(getElement());
bstpshape_m = element_ref_m->getType();
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(element_ref_m->getName(), !Options::asciidump));
bstp_m = dynamic_cast<BeamStripping*>(getElement());
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
const gsl_rng_type* T;
gsl_rng_env_setup();
......@@ -100,9 +100,6 @@ BeamStrippingPhysics::~BeamStrippingPhysics() {
gsl_rng_free(r_m);
}
const std::string BeamStrippingPhysics::getType() const {
return "BeamStrippingPhysics";
}
void BeamStrippingPhysics::apply(PartBunchBase<double, 3>* bunch,
const std::pair<Vector_t, double>& /*boundingSphere*/) {
......@@ -624,7 +621,7 @@ void BeamStrippingPhysics::print(Inform& /*msg*/) {
}
bool BeamStrippingPhysics::stillActive() {
return locPartsInMat_m != 0;
return totalPartsInMat_m != 0;
}
/*
......
......@@ -32,8 +32,8 @@
#include <gsl/gsl_rng.h>
template <class T, unsigned Dim>
class PartBunchBase;
class LossDataSink;
class Inform;
class Cyclotron;
......@@ -55,11 +55,11 @@ public:
virtual void print(Inform& msg);
virtual bool stillActive();
virtual inline double getTime() {return T_m;}
virtual std::string getName() {return element_ref_m->getName();}
virtual size_t getParticlesInMat() {return locPartsInMat_m;}
virtual unsigned getRediffused() {return rediffusedStat_m;}
virtual unsigned int getNumEntered() {return bunchToMatStat_m;}
virtual double getTime();
virtual std::string getName();
virtual size_t getParticlesInMat();
virtual unsigned getRediffused();
virtual unsigned int getNumEntered();
inline void doPhysics(PartBunchBase<double, 3>* bunch);
......@@ -93,7 +93,6 @@ private:
Cyclotron* cycl_m;
BeamStripping* bstp_m;
ElementBase::ElementType bstpshape_m;
gsl_rng* r_m;
......@@ -113,7 +112,7 @@ private:
unsigned bunchToMatStat_m;
unsigned stoppedPartStat_m;
unsigned rediffusedStat_m;
size_t locPartsInMat_m;
size_t totalPartsInMat_m;
static const double csCoefSingle_Hminus[3][9];
static const double csCoefDouble_Hminus[3][9];
......@@ -137,4 +136,34 @@ private:
static double b_m[3][9];
};
inline
double BeamStrippingPhysics::getTime() {
return T_m;
}
inline
std::string BeamStrippingPhysics::getName() {
return (element_ref_m->getName() + "_" + name_m);
}
inline
size_t BeamStrippingPhysics::getParticlesInMat() {
return totalPartsInMat_m;
}
inline
unsigned int BeamStrippingPhysics::getRediffused() {
return rediffusedStat_m;
}
inline
unsigned int BeamStrippingPhysics::getNumEntered() {
return bunchToMatStat_m;
}
inline
const std::string BeamStrippingPhysics::getType() const {
return "BeamStrippingPhysics";
}
#endif //BEAMSTRIPPINGPHYSICS_HH
......@@ -109,6 +109,11 @@ CollimatorPhysics::CollimatorPhysics(const std::string& name,
A3_c(0.0),
A4_c(0.0),
A5_c(0.0),
B1_c(0.0),
B2_c(0.0),
B3_c(0.0),
B4_c(0.0),
B5_c(0.0),
bunchToMatStat_m(0),
stoppedPartStat_m(0),
rediffusedStat_m(0),
......@@ -165,9 +170,6 @@ CollimatorPhysics::~CollimatorPhysics() {
gsl_rng_free(rGen_m);
}
std::string CollimatorPhysics::getName() {
return element_ref_m->getName();
}
/// The material of the collimator
// ------------------------------------------------------------------------
......@@ -215,9 +217,8 @@ void CollimatorPhysics::apply(PartBunchBase<double, 3>* bunch,
dT_m = bunch->getdT();
T_m = bunch->getT();
mass_m = bunch->getM();
charge_m = bunch->getQ();
mass_m = bunch->getM();
charge_m = bunch->getQ();
bool onlyOneLoopOverParticles = ! (allParticleInMat_m);
......@@ -289,7 +290,7 @@ void CollimatorPhysics::computeInteraction(PartBunchBase<double, 3>* bunch) {
// The particle is stopped in the material, set label to -1
locParts_m[i].label = -1.0;
++ stoppedPartStat_m;
lossDs_m->addParticle(R, P, -locParts_m[i].IDincol);
lossDs_m->addParticle(R, P, locParts_m[i].IDincol);
}
}
}
......@@ -305,7 +306,7 @@ void CollimatorPhysics::computeInteraction(PartBunchBase<double, 3>* bunch) {
}
/// Energy Loss: using the Bethe-Bloch equation.
/// In low-energy region use Andersen-Ziegler fitting (only for protons)
/// In low-energy region use Andersen-Ziegler fitting (only for protons and alpha)
/// Energy straggling: For relatively thick absorbers such that the number of collisions
/// is large, the energy loss distribution is shown to be Gaussian in form.
/// See Particle Physics Booklet, chapter 'Passage of particles through matter' or
......@@ -553,7 +554,7 @@ void CollimatorPhysics::copyFromBunch(PartBunchBase<double, 3>* bunch,
// location form the edge of the material. Only use this time step
// for the computation of the interaction with the material, not for
// the integration of the particles. This will ensure that the momenta
// of all particles are reduced by approcimately the same amount in
// of all particles are reduced by approximately the same amount in
// computeEnergyLoss.
double betaz = bunch->P[i](2) / Util::getGamma(bunch->P[i]);
double stepWidth = betaz * Physics::c * bunch->dt[i];
......@@ -603,8 +604,8 @@ void CollimatorPhysics::print(Inform &msg) {
OPALTimer::Timer time;
msg << level2
<< "--- CollimatorPhysics - Name " << element_ref_m->getName()
<< " Material " << material_m << "\n"
<< "--- CollimatorPhysics - Name: " << name_m << "\n"
<< "Material: " << material_m << " - Element: " << element_ref_m->getName() << "\n"
<< "Particle Statistics @ " << time.time() << "\n"
<< std::setw(21) << "entered: " << Util::toStringWithThousandSep(bunchToMatStat_m) << "\n"
<< std::setw(21) << "rediffused: " << Util::toStringWithThousandSep(rediffusedStat_m) << "\n"
......
......@@ -75,7 +75,6 @@ public:
const std::pair<Vector_t, double>& boundingSphere);
virtual const std::string getType() const;
virtual void print(Inform& os);
virtual bool stillActive();
......@@ -84,6 +83,7 @@ public:
virtual size_t getParticlesInMat();
virtual unsigned getRediffused();
virtual unsigned int getNumEntered();
void computeInteraction(PartBunchBase<double, 3>* bunch);
virtual bool computeEnergyLoss(PartBunchBase<double, 3>* bunch,
......@@ -123,8 +123,7 @@ private:
double mass_m; // mass from bunch (eV)
double charge_m; // charge from bunch (elementary charges)
gsl_rng *rGen_m; // random number generator
gsl_rng* rGen_m; // random number generator
std::string material_m; // type of material e.g. aluminum
std::unique_ptr<InsideTester> hitTester_m; // tests whether particles are inside material
ElementBase::ElementType collshape_m; // the type of element (DEGRADER, CCOLLIMATOR or FLEXIBLECOLLIMATOR)
......@@ -193,6 +192,11 @@ double CollimatorPhysics::getTime() {
return T_m;
}
inline
std::string CollimatorPhysics::getName() {
return (element_ref_m->getName() + "_" + name_m);
}
inline
size_t CollimatorPhysics::getParticlesInMat() {
return totalPartsInMat_m;
......
......@@ -51,11 +51,10 @@ public:
Vector_t& P,
const double deltat,
bool includeFluctuations = true) const = 0;
protected:
ElementBase* element_ref_m;
bool allParticleInMat_m; ///< if all particles are in matter stay inside the particle matter interaction
private:
const std::string name_m;
};
......
......@@ -27,9 +27,6 @@
#include "Structure/ParticleMatterInteraction.h"
// Class OpalBeamStripping
// ------------------------------------------------------------------------
OpalBeamStripping::OpalBeamStripping():
OpalElement(SIZE, "BEAMSTRIPPING",
"The \"BEAMSTRIPPING\" element defines a beam stripping interaction"),
......@@ -53,7 +50,7 @@ OpalBeamStripping::OpalBeamStripping():
}
OpalBeamStripping::OpalBeamStripping(const std::string &name, OpalBeamStripping *parent):
OpalBeamStripping::OpalBeamStripping(const std::string& name, OpalBeamStripping* parent):
OpalElement(name, parent),
parmatint_m(NULL) {
setElement(new BeamStrippingRep(name));
......@@ -65,7 +62,7 @@ OpalBeamStripping::~OpalBeamStripping() {
}
OpalBeamStripping *OpalBeamStripping::clone(const std::string &name) {
OpalBeamStripping* OpalBeamStripping::clone(const std::string& name) {
return new OpalBeamStripping(name, this);
}
......@@ -73,8 +70,8 @@ OpalBeamStripping *OpalBeamStripping::clone(const std::string &name) {
void OpalBeamStripping::update() {
OpalElement::update();
BeamStrippingRep *bstp =
dynamic_cast<BeamStrippingRep *>(getElement());
BeamStrippingRep* bstp =
dynamic_cast<BeamStrippingRep*>(getElement());
double pressure = Attributes::getReal(itsAttr[PRESSURE]);
double temperature = Attributes::getReal(itsAttr[TEMPERATURE]);
......@@ -91,7 +88,9 @@ void OpalBeamStripping::update() {
bstp->setResidualGas(gas);
if(itsAttr[PARTICLEMATTERINTERACTION] && parmatint_m == NULL) {
parmatint_m = (ParticleMatterInteraction::find(Attributes::getString(itsAttr[PARTICLEMATTERINTERACTION])))->clone(getOpalName() + std::string("_parmatint"));
const std::string matterDescriptor = Attributes::getString(itsAttr[PARTICLEMATTERINTERACTION]);
ParticleMatterInteraction* orig = ParticleMatterInteraction::find(matterDescriptor);
parmatint_m = orig->clone(matterDescriptor);
parmatint_m->initParticleMatterInteractionHandler(*bstp);
bstp->setParticleMatterInteraction(parmatint_m->handler_m);
}
......
......@@ -25,9 +25,6 @@
#include "Elements/OpalElement.h"
class ParticleMatterInteraction;
// Class OpalBeamStripping
// ------------------------------------------------------------------------
/// The BEAMSTRIPPING element.
class OpalBeamStripping: public OpalElement {
......@@ -50,7 +47,7 @@ public:
virtual ~OpalBeamStripping();
/// Make clone.
virtual OpalBeamStripping *clone(const std::string &name);
virtual OpalBeamStripping* clone(const std::string& name);
/// Update the embedded CLASSIC beam stripping.
virtual void update();
......@@ -58,12 +55,12 @@ public:
private:
// Not implemented.
OpalBeamStripping(const OpalBeamStripping &);
void operator=(const OpalBeamStripping &);
OpalBeamStripping(const OpalBeamStripping&);
void operator=(const OpalBeamStripping&);
// Clone constructor.
OpalBeamStripping(const std::string &name, OpalBeamStripping *parent);
ParticleMatterInteraction *parmatint_m;
OpalBeamStripping(const std::string& name, OpalBeamStripping* parent);
ParticleMatterInteraction* parmatint_m;
};
#endif // OPAL_OpalBeamStripping_HH
\ No newline at end of file
......@@ -22,29 +22,26 @@
#include "Physics/Physics.h"
// Class OpalCCollimator
// ------------------------------------------------------------------------
OpalCCollimator::OpalCCollimator():
OpalElement(SIZE, "CCOLLIMATOR",
"The \"CCOLLIMATOR\" element defines a rectangular-shape cyclotron collimator"),
parmatint_m(NULL) {
itsAttr[XSTART] = Attributes::makeReal
("XSTART", " Start of x coordinate [mm]");
itsAttr[XEND] = Attributes::makeReal
("XEND", " End of x coordinate, [mm]");
itsAttr[XEND] = Attributes::makeReal
("XEND", " End of x coordinate, [mm]");
itsAttr[YSTART] = Attributes::makeReal
("YSTART", "Start of y coordinate, [mm]");
itsAttr[YEND] = Attributes::makeReal
("YEND", "End of y coordinate, [mm]");
itsAttr[YEND] = Attributes::makeReal
("YEND", "End of y coordinate, [mm]");
itsAttr[ZSTART] = Attributes::makeReal
("ZSTART", "Start of vertical coordinate, [mm], default value: -100",-100.0);
itsAttr[ZEND] = Attributes::makeReal
("ZEND", "End of vertical coordinate, [mm], default value: 100", 100.0);
itsAttr[WIDTH] = Attributes::makeReal
("WIDTH", "Width of the collimator [mm]");
itsAttr[OUTFN] = Attributes::makeString
("OUTFN", "Output filename");
("ZSTART", "Start of vertical coordinate, [mm], default value: -100",-100.0);
itsAttr[ZEND] = Attributes::makeReal
("ZEND", "End of vertical coordinate, [mm], default value: 100", 100.0);
itsAttr[WIDTH] = Attributes::makeReal
("WIDTH", "Width of the collimator [mm]");
itsAttr[OUTFN] = Attributes::makeString
("OUTFN", "Output filename");
registerOwnership();
......@@ -52,7 +49,7 @@ OpalCCollimator::OpalCCollimator():
}
OpalCCollimator::OpalCCollimator(const std::string &name, OpalCCollimator *parent):
OpalCCollimator::OpalCCollimator(const std::string& name, OpalCCollimator* parent):
OpalElement(name, parent),
parmatint_m(NULL) {
setElement(new CCollimatorRep(name));
......@@ -60,12 +57,12 @@ OpalCCollimator::OpalCCollimator(const std::string &name, OpalCCollimator *paren
OpalCCollimator::~OpalCCollimator() {
if(parmatint_m)
if (parmatint_m)
delete parmatint_m;
}
OpalCCollimator *OpalCCollimator::clone(const std::string &name) {
OpalCCollimator* OpalCCollimator::clone(const std::string& name) {
return new OpalCCollimator(name, this);
}
......@@ -73,8 +70,9 @@ OpalCCollimator *OpalCCollimator::clone(const std::string &name) {
void OpalCCollimator::update() {
OpalElement::update();
CCollimatorRep *coll =
dynamic_cast<CCollimatorRep *>(getElement());
CCollimatorRep* coll =
dynamic_cast<CCollimatorRep*>(getElement());
const double mm2m = 1e-3;
double xstart = mm2m * Attributes::getReal(itsAttr[XSTART]);
double xend = mm2m * Attributes::getReal(itsAttr[XEND]);
......@@ -90,12 +88,14 @@ void OpalCCollimator::update() {
coll->setDimensions(xstart, xend, ystart, yend, zstart, zend, width);
coll->setOutputFN(Attributes::getString(itsAttr[OUTFN]));
if(itsAttr[PARTICLEMATTERINTERACTION] && parmatint_m == NULL) {
parmatint_m = (ParticleMatterInteraction::find(Attributes::getString(itsAttr[PARTICLEMATTERINTERACTION])))->clone(getOpalName() + std::string("_parmatint"));
if (itsAttr[PARTICLEMATTERINTERACTION] && parmatint_m == NULL) {
const std::string matterDescriptor = Attributes::getString(itsAttr[PARTICLEMATTERINTERACTION]);
ParticleMatterInteraction* orig = ParticleMatterInteraction::find(matterDescriptor);
parmatint_m = orig->clone(matterDescriptor);
parmatint_m->initParticleMatterInteractionHandler(*coll);
coll->setParticleMatterInteraction(parmatint_m->handler_m);
}
// Transmit "unknown" attributes.
OpalElement::updateUnknown(coll);
}
\ No newline at end of file
}
......@@ -45,7 +45,7 @@ public:
virtual ~OpalCCollimator();
/// Make clone.
virtual OpalCCollimator *clone(const std::string &name);
virtual OpalCCollimator* clone(const std::string& name);
/// Update the embedded CLASSIC collimator.
virtual void update();
......@@ -53,12 +53,12 @@ public:
private:
// Not implemented.
OpalCCollimator(const OpalCCollimator &);
void operator=(const OpalCCollimator &);
OpalCCollimator(const OpalCCollimator&);
void operator=(const OpalCCollimator&);
// Clone constructor.
OpalCCollimator(const std::string &name, OpalCCollimator *parent);
ParticleMatterInteraction *parmatint_m;
OpalCCollimator(const std::string& name, OpalCCollimator* parent);
ParticleMatterInteraction* parmatint_m;
};
#endif // OPAL_OpalCCollimator_HH
......@@ -38,7 +38,7 @@ OpalDegrader::OpalDegrader():
}
OpalDegrader::OpalDegrader(const std::string &name, OpalDegrader *parent):
OpalDegrader::OpalDegrader(const std::string& name, OpalDegrader* parent):
OpalElement(name, parent),
parmatint_m(NULL) {
setElement(new DegraderRep(name));
......@@ -46,12 +46,12 @@ OpalDegrader::OpalDegrader(const std::string &name, OpalDegrader *parent):
OpalDegrader::~OpalDegrader() {
if(parmatint_m)
if (parmatint_m)
delete parmatint_m;
}
OpalDegrader *OpalDegrader::clone(const std::string &name) {
OpalDegrader* OpalDegrader::clone(const std::string& name) {
return new OpalDegrader(name, this);
}
......@@ -59,19 +59,22 @@ OpalDegrader *OpalDegrader::clone(const std::string &name) {
void OpalDegrader::update() {
OpalElement::update();
DegraderRep *deg =
dynamic_cast<DegraderRep *>(getElement());
DegraderRep* deg =
dynamic_cast<DegraderRep*>(getElement());
double length = Attributes::getReal(itsAttr[LENGTH]);
deg->setElementLength(length);
deg->setOutputFN(Attributes::getString(itsAttr[OUTFN]));
if(itsAttr[PARTICLEMATTERINTERACTION] && parmatint_m == NULL) {
parmatint_m = (ParticleMatterInteraction::find(Attributes::getString(itsAttr[PARTICLEMATTERINTERACTION])))->clone(getOpalName() + std::string("_parmatint"));
if (itsAttr[PARTICLEMATTERINTERACTION] && parmatint_m == NULL) {
const std::string matterDescriptor = Attributes::getString(itsAttr[PARTICLEMATTERINTERACTION]);
ParticleMatterInteraction* orig = ParticleMatterInteraction::find(matterDescriptor);
parmatint_m = orig->clone(matterDescriptor);
parmatint_m->initParticleMatterInteractionHandler(*deg);
deg->setParticleMatterInteraction(parmatint_m->handler_m);
}
// Transmit "unknown" attributes.
OpalElement::updateUnknown(deg);
}
\ No newline at end of file
}
......@@ -40,7 +40,7 @@ public:
virtual ~OpalDegrader();
/// Make clone.
virtual OpalDegrader *clone(const std::string &name);
virtual OpalDegrader* clone(const std::string& name);
/// Update the embedded CLASSIC collimator.
virtual void update();
......@@ -48,13 +48,13 @@ public:
private:
// Not implemented.
OpalDegrader(const OpalDegrader &);
void operator=(const OpalDegrader &);
OpalDegrader(const OpalDegrader&);
void operator=(const OpalDegrader&);
// Clone constructor.
OpalDegrader(const std::string &name, OpalDegrader *parent);
OpalDegrader(const std::string& name, OpalDegrader* parent);
ParticleMatterInteraction *parmatint_m;
ParticleMatterInteraction* parmatint_m;
};
#endif // OPAL_OpalDegrader_HH
#endif // OPAL_OpalDegrader_HH
\ No newline at end of file
// ------------------------------------------------------------------------
// $RCSfile: OpalDrift.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: OpalDrift
// Class OpalDrift
// The class of OPAL drift spaces.
//
// ------------------------------------------------------------------------
// Copyright (c) 200x - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// $Date: 2000/03/27 09:33:39 $
// $Author: Andreas Adelmann $
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
// ------------------------------------------------------------------------
#include "Elements/OpalDrift.h"
#include "Structure/BoundaryGeometry.h"
#include "Attributes/Attributes.h"
......@@ -23,8 +22,6 @@
#include "Structure/OpalWake.h"
#include "Structure/ParticleMatterInteraction.h"
// Class OpalDrift
// ------------------------------------------------------------------------
OpalDrift::OpalDrift():
OpalElement(SIZE, "DRIFT",
......@@ -41,9 +38,8 @@ OpalDrift::OpalDrift():
itsAttr[GEOMETRY] = Attributes::makeString
("GEOMETRY", "BoundaryGeometry for Drifts");
itsAttr[NSLICES] = Attributes::makeReal
("NSLICES",
"The number of slices/ steps for this element in Map Tracking", 1);
itsAttr[NSLICES] = Attributes::makeReal
("NSLICES", "The number of slices/ steps for this element in Map Tracking", 1);
registerOwnership();
......@@ -51,7 +47,7 @@ OpalDrift::OpalDrift():
}
OpalDrift::OpalDrift(const std::string &name, OpalDrift *parent):
OpalDrift::OpalDrift(const std::string& name, OpalDrift* parent):
OpalElement(name, parent),
owk_m(NULL),
parmatint_m(NULL),
......@@ -61,16 +57,16 @@ OpalDrift::OpalDrift(const std::string &name, OpalDrift *parent):
OpalDrift::~OpalDrift() {
if(owk_m)
if (owk_m)
delete owk_m;
if(parmatint_m)
if (parmatint_m)
delete parmatint_m;
if(obgeo_m)
if (obgeo_m)
delete obgeo_m;
}
OpalDrift *OpalDrift::clone(const std::string &name) {
OpalDrift* OpalDrift::clone(const std::string& name) {
return new OpalDrift(name, this);
}
......@@ -83,23 +79,28 @@ bool OpalDrift::isDrift() const {
void OpalDrift::update() {
OpalElement::update();
DriftRep *drf = static_cast<DriftRep *>(getElement());
DriftRep* drf = static_cast<DriftRep*>(getElement());
drf->setElementLength(Attributes::getReal(itsAttr[LENGTH]));
drf->setNSlices(Attributes::getReal(itsAttr[NSLICES]));
if(itsAttr[WAKEF] && owk_m == NULL) {
if (itsAttr[WAKEF] && owk_m == NULL) {
owk_m = (OpalWake::find(Attributes::getString(itsAttr[WAKEF])))->clone(getOpalName() + std::string("_wake"));