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