KaliVeda
Toolkit for HIC analysis
KVElasticCountRates.cpp
1 /*
2 $Id: KVElasticCountRates.cpp,v 1.9 2007/04/04 10:39:17 ebonnet Exp $
3 $Revision: 1.9 $
4 $Date: 2007/04/04 10:39:17 $
5 */
6 
7 //Created by KVClassFactory on Fri Nov 20 2015
8 //Author: John Frankland
9 
10 #include "KVElasticCountRates.h"
11 #include "KVMultiDetArray.h"
12 #include "TH1F.h"
13 #include "KVDetector.h"
14 #include "KVGroup.h"
15 #include "KVTarget.h"
16 #include "KV2Body.h"
17 #include "KVNucleus.h"
18 #include "TObjArray.h"
19 #include "KVDatime.h"
20 
21 #include <TH2F.h>
22 #include <vector>
23 #include <algorithm>
24 
25 using namespace std;
26 
29 
30 
31 
34  const TVector3& beam_dir):
35  fAngularRange(theta_min, theta_max, phi_min, phi_max),
36  fBeamDirection{beam_dir}
37 {
38  Fatal("KVElasticCountRates", "This class must be reimplemented using KVDetectionSimulator for detection of elastically scattered particles.");
39  //Default constructor
40  fDepth = fTheta = 0;
41  fBinE = 500;
42  fEnergy = 0;
43  fKinematics = 0;
44  fTarget = 0;
45  fIntLayer = 0;
46  fMultiLayer = kFALSE;
47  fVolume = (theta_max - theta_min) * (phi_max - phi_min) * TMath::DegToRad() * TMath::DegToRad();
48  //fVolume = (cos(theta_min*TMath::DegToRad())-cos(theta_max*TMath::DegToRad()))*(phi_max-phi_min)*TMath::DegToRad();
49  fExx = 0.;
50  if (gMultiDetArray) {
51  gMultiDetArray->SetSimMode(kTRUE);
53  gMultiDetArray->InitializeIDTelescopes();
54  }
55  else {
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)");
58  }
59 }
60 
61 
62 
63 
66 
67 KVElasticCountRates::~KVElasticCountRates()
68 {
69  //Destructor
70  if (fDepth)
71  delete fDepth;
72  if (fTheta)
73  delete fTheta;
74  if (fKinematics)
75  delete fKinematics;
76 }
77 
78 
79 
80 
83 
85 {
86  //Set detector parameters, target, etc. for run
87  gMultiDetArray->SetParameters(run);
88  gMultiDetArray->InitializeIDTelescopes();
89  fTarget = gMultiDetArray->GetTarget();
90  fAtomicDensity = fTarget->GetAtomsPerCM2() * 1.e-24; //in barn^-1 units
91  fTarget->SetRandomized();
92  fIntLayer = 0;
93  fMultiLayer = (fTarget->NumberOfLayers() > 1);
94 }
95 
96 
97 
100 
102 {
103  // Set projectile and beam energy [MeV/nucleon] & direction
104  fProj = KVNucleus(nuc, e_sur_a);
105  fEnergy = fProj.GetE();
106  fProj.SetMomentum(fEnergy, fBeamDirection);
107 }
108 
109 
110 
111 
117 
119 {
120  //For multilayer targets, use this method to choose in which
121  //layer the scattering will take place.
122  //If name="", reset any previous choice so that scattering
123  //can take place in any layer
124  if (!fTarget) {
125  cout <<
126  "<KVElasticCountRates::SetTargetScatteringLayer> : No target set. Set run first."
127  << endl;
128  return;
129  }
130  fIntLayer = fTarget->GetLayerIndex(name);
131  if (fIntLayer)
132  fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
133 }
134 
135 
136 
137 
142 
144 {
145  //Set the number of bins of the GetEnergy() histogram
146  //Default value is 500; this function has to be called before
147  //using CalculateScattering.
148  fBinE = nbins;
149 }
150 
151 
152 
158 
160 {
161  //Perform scattering 'N' times for current values
162  //of particle Z, A and energy, target excited state, and detector.
163  //Print out for each hit detector the mean energy loss & counting rate
164  //for a beam intensity of 10**7 pps
165 
166  fNtirages = N;
167 
168  fHistos.Delete();
169  fRates.clear();
170 
171  if (!fProj.IsDefined()) {
172  cout <<
173  "<KVElasticCountRates::CalculateScattering> : Set projectile properties first"
174  << endl;
175  return;
176  }
177  if (!fEnergy) {
178  cout <<
179  "<KVElasticCountRates::CalculateScattering> : Set projectile energy first"
180  << endl;
181  return;
182  }
183  if (!fTarget) {
184  cout <<
185  "<KVElasticCountRates::CalculateScattering> : No target set. Set run first."
186  << endl;
187  return;
188  }
189 
190  /* -------------------------------------------------------------------------------------------------------------------------- */
191 
192  //set up histograms
193 
194  /* -------------------------------------------------------------------------------------------------------------------------- */
195 
196  // -- histograms for debugging: depth in target and angular distribution
197  if (fDepth)
198  delete fDepth;
199  if (fTheta)
200  delete fTheta;
201  fDepth =
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());
207 
208  /* -------------------------------------------------------------------------------------------------------------------------- */
209 
210  //set up kinematics
211  if (!fKinematics)
212  fKinematics = new KV2Body;
213 
214  /* -------------------------------------------------------------------------------------------------------------------------- */
215 
216  //set random interaction point for scattering
217  if (fIntLayer) {
218  fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
219  }
220  else {
221  fTarget->GetInteractionPoint(&fProj);
222  }
223 
224  /* -------------------------------------------------------------------------------------------------------------------------- */
225 
226  //get target nucleus properties from scattering layer
227  TVector3 IP = fTarget->GetInteractionPoint();
228  KVMaterial* targ_mat = fTarget->GetLayer(IP);
229  KVNucleus t;
230  t.SetZandA((Int_t) targ_mat->GetZ(), (Int_t) targ_mat->GetMass());
231  fKinematics->SetTarget(t);
232 
233  /* -------------------------------------------------------------------------------------------------------------------------- */
234 
235  //set excited state of target nucleus - in other words, dissipated energy for
236  //reaction due to excitation of target
237  fKinematics->SetEDiss(fExx);
238 
239  /* -------------------------------------------------------------------------------------------------------------------------- */
240 
241  for (int i = 0; i < N; i++) {
242  //calculate slowing of incoming projectile
243  fTarget->SetIncoming();
244  fTarget->DetectParticle(&fProj);
245  fKinematics->SetProjectile(fProj);
246  fKinematics->SetOutgoing(fProj);
247  fKinematics->CalculateKinematics();
248  //set random direction of outgoing projectile
249  fAngularRange.GetRandomAngles(theta, phi, "random");
250  fProj.SetTheta(theta);
251  fProj.SetPhi(phi);
252  // get Rutherford Xsec for projectile scattered to theta degrees in lab
253  // taking into account the initial direction of the beam (if beam direction = (0,0,1)
254  // this is just the same as theta)
255  double angle_with_beam = TMath::RadToDeg() * fProj.Angle(fBeamDirection);
256  xsec = TMath::Abs(fKinematics->GetXSecRuthLab(angle_with_beam));
257  //set energy of scattered nucleus
258  //WARNING: for inverse kinematics reactions, their are two energies for
259  //each angle below the maximum scattering angle.
260  //We only use the highest energy corresponding to the most forward CM angle.
261  Double_t e1, e2;
262  fKinematics->GetELab(3, angle_with_beam, 3, e1, e2);
263  fProj.SetEnergy(TMath::Max(e1, e2));
264  fTheta->Fill(theta, xsec);
265  //slowing of outgoing projectile in target
266  fTarget->SetOutgoing();
267  fTarget->DetectParticle(&fProj);
268  //now detect particle in array
269  KVNameValueList* detectors = nullptr;//gMultiDetArray->DetectParticle(&fProj);
270  //fill histograms
271  fDepth->Fill(IP.z());
272  FillHistograms(detectors);
273  fProj.GetParameters()->Clear();
274  fProj.SetEnergy(fEnergy);
275  fProj.SetMomentum(fEnergy, fBeamDirection);
276  //set random interaction point for scattering
277  if (fIntLayer) {
278  fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
279  }
280  else {
281  fTarget->GetInteractionPoint(&fProj);
282  //if target is multilayer and the interaction layer is not fixed,
283  //the layer & hence the target nucleus may change
284  if (fMultiLayer) {
285  targ_mat = fTarget->GetLayer(fTarget->GetInteractionPoint());
286  KVNucleus t;
287  t.SetZandA((Int_t) targ_mat->GetZ(), (Int_t) targ_mat->GetMass());
288  fKinematics->SetTarget(t);
289  }
290  }
291  IP = fTarget->GetInteractionPoint();
292  }
293 
294  PrintResults();
295 }
296 
297 
298 
299 
300 
301 
307 
309 {
310  // parse the list dets
311  // fill histograms with energy loss for all detectors
312  // clear the detector energy losses
313  // delete the list
314 
315  if (!dets) return;
316 
317  Int_t ndets = dets->GetNpar();
318  for (int i = 0; i < ndets; i++) {
319  TString detname = dets->GetNameAt(i);
320  KVDetector* det = gMultiDetArray->GetDetector(detname);
321  if (!det) continue;
322  TH1F* histo = (TH1F*)fHistos.FindObject(detname);
323  if (!histo) {
324  histo = new TH1F(detname, Form("Eloss in %s", detname.Data()), fBinE, 0, 0);
325  fHistos.Add(histo);
326  }
327  double de = dets->GetDoubleValue(i);
328  histo->Fill(de, xsec * sin(theta * TMath::DegToRad()));
329  histo = (TH1F*)fHistos.FindObject(detname + "_dW");
330  if (!histo) {
331  histo = new TH1F(detname + "_dW", Form("Solid angle of %s", detname.Data()), fBinE, 0, 0);
332  fHistos.Add(histo);
333  }
334  histo->Fill(de, sin(theta * TMath::DegToRad()));
335  TH2F* histo2 = (TH2F*)fHistos.FindObject(detname + "_map");
336  if (!histo2) {
337  histo2 = new TH2F(detname + "_map", Form("Map of %s", detname.Data()), 100, 0, 0, 100, 0, 0);
338  fHistos.Add(histo2);
339  }
340  histo2->Fill(theta, phi, xsec);
341  fGlobalMap->Fill(theta * sin(phi * TMath::DegToRad()), theta * cos(phi * TMath::DegToRad()), xsec);
342  det->Clear();
343  }
344  delete dets;
345 }
346 
347 
351 
356 struct count_rate {
358  double counts;
359  double energy;
360  double theta;
361  double phi;
362  double fluence;
363  double dissipation;
364  double intXsec;
365  count_rate(TString n, double c, double e, double t, double p, double f, double d, double i)
366  : detector(n), counts(c), energy(e), theta(t), phi(p), fluence(f), dissipation(d), intXsec(i) {}
367  void print()
368  {
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);
371  }
373  {
374  if (nl.HasParameter("DetList")) {
375  KVString tmp = nl.GetStringValue("DetList");
376  tmp += ",";
377  tmp += detector;
378  nl.SetValue("DetList", tmp.Data());
379  }
380  else {
381  nl.SetValue("DetList", detector.Data());
382  }
383  nl.SetValue(Form("%s.Counts(/sec)", detector.Data()), Form("%8.2f", counts));
384 
385  nl.SetValue(Form("%s.Counts(/sec)", detector.Data()), Form("%8.2f", counts));
386  nl.SetValue(Form("%s.Emean(MeV)", detector.Data()), Form("%7.1f", energy));
387  nl.SetValue(Form("%s.TotXsec(barn)", detector.Data()), Form("%9.3E", intXsec));
388  nl.SetValue(Form("%s.Fluence(/sec./cm**2)", detector.Data()), Form("%9.3E", fluence));
389  nl.SetValue(Form("%s.dissip(MeV/sec./cm**2)", detector.Data()), Form("%9.3E", dissipation));
390  }
391 };
392 
393 
395 
396 bool compare_count_rates(count_rate a, count_rate b)
397 {
398  if (TMath::Abs(a.theta - b.theta) < 0.5) return a.phi < b.phi;
399  return a.theta < b.theta;
400 }
401 
402 
403 
406 
408 {
409  // Print mean energy deposit & counting rate for given beam intensity in particles per second
410 
411  TIter it(&fHistos);
412  TH1F* h;
413  fRates.clear();
414 
415  std::vector<count_rate> count_rates;
416 
417  while ((h = (TH1F*)it())) {
418  TString name = h->GetName();
419  if (!name.EndsWith("_dW") && !name.EndsWith("_map")) {
420  TH2F* map = (TH2F*)fHistos.FindObject(name + "_map");
421  // integrated cross-section
422  double intXsec = h->Integral() * fVolume / fNtirages;
423  // counting rate
424  double rate = fAtomicDensity * beam_intensity * intXsec;
425  // mean energy
426  double emean = h->GetMean();
427  KVDetector* det = gMultiDetArray->GetDetector(name);
428  double fluence = rate / det->GetEntranceWindowSurfaceArea();
429  double dissipation = emean * rate / det->GetEntranceWindowSurfaceArea();
430  count_rates.push_back(
431  count_rate(name, rate, emean, map->GetMean(), map->GetMean(2), fluence, dissipation, intXsec)
432  );
433  fRates[name.Data()] = KVElasticCountRate(rate, emean, intXsec, fluence, dissipation);
434  }
435  }
436  std::sort(count_rates.begin(), count_rates.end(), compare_count_rates);
437 
438  for (std::vector<count_rate>::iterator it = count_rates.begin(); it != count_rates.end(); ++it) {
439  it->print();
440  }
441 }
442 
443 
444 
447 
449 {
450  // Print mean energy deposit & counting rate for given beam intensity in particles per second
451  KVNameValueList nl;
452  nl.SetName("Generated by KVElasticCountRates::PutResultsInList method");
453  KVDatime dt;
454  nl.SetTitle(dt.AsSQLString());
455  TIter it(&fHistos);
456  TH1F* h;
457  fRates.clear();
458  nl.SetValue("Intensity(pps)", beam_intensity);
459 
460  std::vector<count_rate> count_rates;
461 
462  while ((h = (TH1F*)it())) {
463  TString name = h->GetName();
464  if (!name.EndsWith("_dW") && !name.EndsWith("_map")) {
465  TH2F* map = (TH2F*)fHistos.FindObject(name + "_map");
466  // integrated cross-section
467  double intXsec = h->Integral() * fVolume / fNtirages;
468  // counting rate
469  double rate = fAtomicDensity * beam_intensity * intXsec;
470  // mean energy
471  double emean = h->GetMean();
472  KVDetector* det = gMultiDetArray->GetDetector(name);
473  double fluence = rate / det->GetEntranceWindowSurfaceArea();
474  double dissipation = emean * rate / det->GetEntranceWindowSurfaceArea();
475  count_rates.push_back(
476  count_rate(name, rate, emean, map->GetMean(), map->GetMean(2), fluence, dissipation, intXsec)
477  );
478  fRates[name.Data()] = KVElasticCountRate(rate, emean, intXsec, fluence, dissipation);
479  }
480  }
481  std::sort(count_rates.begin(), count_rates.end(), compare_count_rates);
482 
483  for (std::vector<count_rate>::iterator it = count_rates.begin(); it != count_rates.end(); ++it) {
484  it->PutInList(nl);
485  }
486  return nl;
487 }
488 
489 
490 //__________________________________________________________________//
491 
int Int_t
#define d(i)
#define f(i)
#define c(i)
#define e(i)
char Char_t
constexpr Bool_t kFALSE
double Double_t
#define N
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 name[80]
char * Form(const char *fmt,...)
Relativistic binary kinematics calculator.
Definition: KV2Body.h:166
Extension of TDatime to handle various useful date formats.
Definition: KVDatime.h:33
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.
Definition: KVMaterial.h:89
Double_t GetZ() const
Definition: KVMaterial.cpp:390
Double_t GetMass() const
Definition: KVMaterial.cpp:302
KVTarget * GetTarget()
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.
Definition: KVNucleus.h:123
void SetZandA(Int_t z, Int_t a)
Set atomic number and mass number.
Definition: KVNucleus.cpp:713
Extension of ROOT TString class which allows backwards compatibility with ROOT v3....
Definition: KVString.h:73
Double_t GetAtomsPerCM2() const
virtual UInt_t GetUnits() const;
Definition: KVTarget.cpp:927
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
Double_t z() const
RVec< PromoteType< T > > cos(const RVec< T > &v)
RVec< PromoteType< T > > sin(const RVec< T > &v)
const Int_t n
TH1 * h
void Fatal(const char *location, const char *fmt,...)
void Warning(const char *location, const char *fmt,...)
constexpr Double_t DegToRad()
Double_t Abs(Double_t d)
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)
TArc a
ClassImp(TPyArg)