10 #include "KVElasticCountRates.h"
11 #include "KVMultiDetArray.h"
13 #include "KVDetector.h"
17 #include "KVNucleus.h"
20 #include "KVDetectionSimulator.h"
37 fAngularRange(theta_min, theta_max, phi_min, phi_max),
38 fBeamDirection{beam_dir}
57 Warning(
"KVElasticCountRates",
"gMultiDetArray does not refer to a valid multidetector array");
58 printf(
"Define it before using this class, and put it in simulation mode : gMultiDetArray->SetSimMode(kTRUE)");
68 KVElasticCountRates::~KVElasticCountRates()
92 fTarget->SetRandomized();
94 fMultiLayer = (fTarget->NumberOfLayers() > 1);
106 fEnergy = fProj.GetE();
107 fProj.SetMomentum(fEnergy, fBeamDirection);
127 "<KVElasticCountRates::SetTargetScatteringLayer> : No target set. Set run first."
131 fIntLayer = fTarget->GetLayerIndex(
name);
133 fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
172 if (!fProj.IsDefined()) {
174 "<KVElasticCountRates::CalculateScattering> : Set projectile properties first"
180 "<KVElasticCountRates::CalculateScattering> : Set projectile energy first"
186 "<KVElasticCountRates::CalculateScattering> : No target set. Set run first."
203 new TH1F(
"hDepth",
"Depth (mg/cm2)", 500, 0., fTarget->GetTotalEffectiveThickness());
204 fTheta =
new TH1F(
"hTheta",
"Theta (deg.)", 500, 0., 0.);
205 fGlobalMap =
new TH2F(
"hGlobalMap",
"X-section in x(vert)-y(hori) plane",
206 500, -fAngularRange.GetThetaMax(), fAngularRange.GetThetaMax(),
207 500, -fAngularRange.GetThetaMax(), fAngularRange.GetThetaMax());
219 fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
222 fTarget->GetInteractionPoint(&fProj);
228 TVector3 IP = fTarget->GetInteractionPoint();
232 fKinematics->SetTarget(t);
238 fKinematics->SetEDiss(fExx);
244 for (
int i = 0; i <
N; i++) {
249 fTarget->SetIncoming();
250 fTarget->DetectParticle(&fProj);
251 fKinematics->SetProjectile(fProj);
252 fKinematics->SetOutgoing(fProj);
253 fKinematics->CalculateKinematics();
255 fAngularRange.GetRandomAngles(theta, phi,
"random");
256 fProj.SetTheta(theta);
261 double angle_with_beam =
TMath::RadToDeg() * fProj.Angle(fBeamDirection);
262 xsec =
TMath::Abs(fKinematics->GetXSecRuthLab(angle_with_beam));
268 fKinematics->GetELab(3, angle_with_beam, 3, e1, e2);
270 fTheta->Fill(theta, xsec);
272 fTarget->SetOutgoing();
273 fTarget->DetectParticle(&fProj);
277 fDepth->Fill(IP.
z());
278 FillHistograms(&detectors);
279 fProj.GetParameters()->Clear();
280 fProj.SetEnergy(fEnergy);
281 fProj.SetMomentum(fEnergy, fBeamDirection);
284 fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
287 fTarget->GetInteractionPoint(&fProj);
291 targ_mat = fTarget->GetLayer(fTarget->GetInteractionPoint());
294 fKinematics->SetTarget(t);
297 IP = fTarget->GetInteractionPoint();
324 for (
int i = 0; i < ndets; i++) {
326 KVDetector* det = gMultiDetArray->
GetDetector(detname);
330 histo =
new TH1F(detname,
Form(
"Eloss in %s", detname.
Data()), fBinE, 0, 0);
337 histo =
new TH1F(detname +
"_dW",
Form(
"Solid angle of %s", detname.
Data()), fBinE, 0, 0);
343 histo2 =
new TH2F(detname +
"_map",
Form(
"Map of %s", detname.
Data()), 100, 0, 0, 100, 0, 0);
346 histo2->
Fill(theta, phi, xsec);
371 : detector(
n), counts(
c), energy(
e), theta(t), phi(
p), fluence(
f), dissipation(
d), intXsec(i) {}
374 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",
375 detector.
Data(), counts, energy, intXsec, fluence, dissipation);
404 return a.theta <
b.theta;
420 std::vector<count_rate> count_rates;
422 while ((
h = (
TH1F*)it())) {
424 if (!
name.EndsWith(
"_dW") && !
name.EndsWith(
"_map")) {
427 double intXsec =
h->
Integral() * fVolume / fNtirages;
429 double rate = fAtomicDensity * beam_intensity * intXsec;
433 double fluence = rate / det->GetEntranceWindowSurfaceArea();
434 double dissipation = emean * rate / det->GetEntranceWindowSurfaceArea();
435 count_rates.push_back(
441 std::sort(count_rates.begin(), count_rates.end(), compare_count_rates);
443 for (std::vector<count_rate>::iterator it = count_rates.begin(); it != count_rates.end(); ++it) {
457 nl.
SetName(
"Generated by KVElasticCountRates::PutResultsInList method");
463 nl.
SetValue(
"Intensity(pps)", beam_intensity);
465 std::vector<count_rate> count_rates;
467 while ((
h = (
TH1F*)it())) {
469 if (!
name.EndsWith(
"_dW") && !
name.EndsWith(
"_map")) {
472 double intXsec =
h->
Integral() * fVolume / fNtirages;
474 double rate = fAtomicDensity * beam_intensity * intXsec;
478 double fluence = rate / det->GetEntranceWindowSurfaceArea();
479 double dissipation = emean * rate / det->GetEntranceWindowSurfaceArea();
480 count_rates.push_back(
486 std::sort(count_rates.begin(), count_rates.end(), compare_count_rates);
488 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.
Simulate detection of events in a detector array.
KVNameValueList PropagateParticle(KVNucleus *)
Calculate elastic scattering count rates in multidetector arrays.
void FillHistograms(const KVNameValueList *)
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 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)