10 #include "KVElasticCountRates.h"
11 #include "KVMultiDetArray.h"
13 #include "KVDetector.h"
17 #include "KVNucleus.h"
35 fAngularRange(theta_min, theta_max, phi_min, phi_max),
36 fBeamDirection{beam_dir}
38 Fatal(
"KVElasticCountRates",
"This class must be reimplemented using KVDetectionSimulator for detection of elastically scattered particles.");
56 Warning(
"KVElasticCountRates",
"gMultiDetArray does not refer to a valid multidetector array");
57 printf(
"Define it before using this class, and put it in simulation mode : gMultiDetArray->SetSimMode(kTRUE)");
67 KVElasticCountRates::~KVElasticCountRates()
91 fTarget->SetRandomized();
93 fMultiLayer = (fTarget->NumberOfLayers() > 1);
105 fEnergy = fProj.GetE();
106 fProj.SetMomentum(fEnergy, fBeamDirection);
126 "<KVElasticCountRates::SetTargetScatteringLayer> : No target set. Set run first."
130 fIntLayer = fTarget->GetLayerIndex(
name);
132 fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
171 if (!fProj.IsDefined()) {
173 "<KVElasticCountRates::CalculateScattering> : Set projectile properties first"
179 "<KVElasticCountRates::CalculateScattering> : Set projectile energy first"
185 "<KVElasticCountRates::CalculateScattering> : No target set. Set run first."
202 new TH1F(
"hDepth",
"Depth (mg/cm2)", 500, 0., fTarget->GetTotalEffectiveThickness());
203 fTheta =
new TH1F(
"hTheta",
"Theta (deg.)", 500, 0., 0.);
204 fGlobalMap =
new TH2F(
"hGlobalMap",
"X-section in x(vert)-y(hori) plane",
205 500, -fAngularRange.GetThetaMax(), fAngularRange.GetThetaMax(),
206 500, -fAngularRange.GetThetaMax(), fAngularRange.GetThetaMax());
218 fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
221 fTarget->GetInteractionPoint(&fProj);
227 TVector3 IP = fTarget->GetInteractionPoint();
231 fKinematics->SetTarget(t);
237 fKinematics->SetEDiss(fExx);
241 for (
int i = 0; i <
N; i++) {
243 fTarget->SetIncoming();
244 fTarget->DetectParticle(&fProj);
245 fKinematics->SetProjectile(fProj);
246 fKinematics->SetOutgoing(fProj);
247 fKinematics->CalculateKinematics();
249 fAngularRange.GetRandomAngles(theta, phi,
"random");
250 fProj.SetTheta(theta);
255 double angle_with_beam =
TMath::RadToDeg() * fProj.Angle(fBeamDirection);
256 xsec =
TMath::Abs(fKinematics->GetXSecRuthLab(angle_with_beam));
262 fKinematics->GetELab(3, angle_with_beam, 3, e1, e2);
264 fTheta->Fill(theta, xsec);
266 fTarget->SetOutgoing();
267 fTarget->DetectParticle(&fProj);
271 fDepth->Fill(IP.
z());
272 FillHistograms(detectors);
273 fProj.GetParameters()->Clear();
274 fProj.SetEnergy(fEnergy);
275 fProj.SetMomentum(fEnergy, fBeamDirection);
278 fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
281 fTarget->GetInteractionPoint(&fProj);
285 targ_mat = fTarget->GetLayer(fTarget->GetInteractionPoint());
288 fKinematics->SetTarget(t);
291 IP = fTarget->GetInteractionPoint();
318 for (
int i = 0; i < ndets; i++) {
320 KVDetector* det = gMultiDetArray->
GetDetector(detname);
324 histo =
new TH1F(detname,
Form(
"Eloss in %s", detname.
Data()), fBinE, 0, 0);
331 histo =
new TH1F(detname +
"_dW",
Form(
"Solid angle of %s", detname.
Data()), fBinE, 0, 0);
337 histo2 =
new TH2F(detname +
"_map",
Form(
"Map of %s", detname.
Data()), 100, 0, 0, 100, 0, 0);
340 histo2->
Fill(theta, phi, xsec);
366 : detector(
n), counts(
c), energy(
e), theta(t), phi(
p), fluence(
f), dissipation(
d), intXsec(i) {}
369 printf(
"%s \t: N=%8.2f/sec. \t <E>=%7.1f MeV \t Tot.Xsec=%9.3E barn \t fluence=%9.3E/sec./cm**2 \t dissip.=%9.3E MeV/sec./cm**2\n",
370 detector.
Data(), counts, energy, intXsec, fluence, dissipation);
399 return a.theta <
b.theta;
415 std::vector<count_rate> count_rates;
417 while ((
h = (
TH1F*)it())) {
419 if (!
name.EndsWith(
"_dW") && !
name.EndsWith(
"_map")) {
422 double intXsec =
h->
Integral() * fVolume / fNtirages;
424 double rate = fAtomicDensity * beam_intensity * intXsec;
428 double fluence = rate / det->GetEntranceWindowSurfaceArea();
429 double dissipation = emean * rate / det->GetEntranceWindowSurfaceArea();
430 count_rates.push_back(
436 std::sort(count_rates.begin(), count_rates.end(), compare_count_rates);
438 for (std::vector<count_rate>::iterator it = count_rates.begin(); it != count_rates.end(); ++it) {
452 nl.
SetName(
"Generated by KVElasticCountRates::PutResultsInList method");
458 nl.
SetValue(
"Intensity(pps)", beam_intensity);
460 std::vector<count_rate> count_rates;
462 while ((
h = (
TH1F*)it())) {
464 if (!
name.EndsWith(
"_dW") && !
name.EndsWith(
"_map")) {
467 double intXsec =
h->
Integral() * fVolume / fNtirages;
469 double rate = fAtomicDensity * beam_intensity * intXsec;
473 double fluence = rate / det->GetEntranceWindowSurfaceArea();
474 double dissipation = emean * rate / det->GetEntranceWindowSurfaceArea();
475 count_rates.push_back(
481 std::sort(count_rates.begin(), count_rates.end(), compare_count_rates);
483 for (std::vector<count_rate>::iterator it = count_rates.begin(); it != count_rates.end(); ++it) {
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
char * Form(const char *fmt,...)
Relativistic binary kinematics calculator.
Extension of TDatime to handle various useful date formats.
Calculate elastic scattering count rates in multidetector arrays.
void PrintResults(Double_t beam_intensity=1.e+07)
Print mean energy deposit & counting rate for given beam intensity in particles per second.
KVNameValueList PutResultsInList(Double_t beam_intensity=1.e+07)
Print mean energy deposit & counting rate for given beam intensity in particles per second.
void FillHistograms(KVNameValueList *)
void SetProjectile(const Char_t *nuc, Double_t e_sur_a)
Set projectile and beam energy [MeV/nucleon] & direction.
void SetEbinning(Int_t nbins=500)
void CalculateScattering(Int_t N=10000)
void SetRun(Int_t run)
Set detector parameters, target, etc. for run.
void SetTargetScatteringLayer(const Char_t *name)
virtual KVDetector * GetDetector(const Char_t *name) const
Return detector in this structure with given name.
Description of physical materials used to construct detectors & targets; interface to range tables.
virtual void SetParameters(UInt_t n, Bool_t physics_parameters_only=kFALSE)
void SetFilterType(Int_t t)
virtual void InitializeIDTelescopes()
virtual void SetSimMode(Bool_t on=kTRUE)
Handles lists of named parameters with different types, a list of KVNamedParameter objects.
Double_t GetDoubleValue(const Char_t *name) const
void SetValue(const Char_t *name, value_type value)
const Char_t * GetNameAt(Int_t idx) const
Int_t GetNpar() const
return the number of stored parameters
const Char_t * GetStringValue(const Char_t *name) const
Bool_t HasParameter(const Char_t *name) const
Description of properties and kinematics of atomic nuclei.
void SetZandA(Int_t z, Int_t a)
Set atomic number and mass number.
Extension of ROOT TString class which allows backwards compatibility with ROOT v3....
Double_t GetAtomsPerCM2() const
virtual UInt_t GetUnits() const;
const char * AsSQLString() const
virtual Double_t GetMean(Int_t axis=1) const
TObject * FindObject(const char *name) const override
virtual Int_t Fill(const char *name, Double_t w)
virtual Double_t Integral(Int_t binx1, Int_t binx2, Option_t *option="") const
virtual Int_t Fill(const char *namex, const char *namey, Double_t w)
virtual void SetTitle(const char *title="")
const char * GetName() const override
virtual void SetName(const char *name)
const char * Data() const
RVec< PromoteType< T > > cos(const RVec< T > &v)
RVec< PromoteType< T > > sin(const RVec< T > &v)
constexpr Double_t DegToRad()
Double_t Max(Double_t a, Double_t b)
constexpr Double_t RadToDeg()
Utility class used by KVElasticCountRates to store results.
Utility class used by KVElasticCountRates.
count_rate(TString n, double c, double e, double t, double p, double f, double d, double i)
void PutInList(KVNameValueList &nl)