KaliVeda
Toolkit for HIC analysis
KVElasticScatter.cpp
1 /*
2 $Id: KVElasticScatter.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 Thu Oct 19 22:31:02 2006
8 //Author: John Frankland
9 
10 #include "KVElasticScatter.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 "TObjArray.h"
18 
19 using namespace std;
20 
22 
23 
24 
27 KVElasticScatter::KVElasticScatter(): fBeamDirection(0, 0, 1)
28 {
29  //Default constructor
30  fDepth = 0;
31  fTheta = 0;
32  fBinE = 500;
33  fEnergy = 0;
34  fKinematics = 0;
35  fTarget = 0;
36  fIntLayer = fNDets = 0;
37  fDetector = 0;
38  fMultiLayer = kFALSE;
39  fExx = 0.;
40  fHistos = 0;
41  if (gMultiDetArray) {
42  gMultiDetArray->SetSimMode(kTRUE);
43  }
44  else {
45  Warning("KVElasticScatter", "gMultiDetArray does not refer to a valid multidetector array");
46  printf("Define it before using this class, and put it in simulation mode : gMultiDetArray->SetSimMode(kTRUE)");
47  }
48 }
49 
50 
51 
52 
55 
56 KVElasticScatter::~KVElasticScatter()
57 {
58  //Destructor
59  if (fDepth)
60  delete fDepth;
61  if (fTheta)
62  delete fTheta;
63  if (fKinematics)
64  delete fKinematics;
65  if (fHistos)
66  delete fHistos;
67  gMultiDetArray->SetSimMode(kFALSE);
68 }
69 
70 
71 
72 
75 
77 {
78  //Set detector parameters, target, etc. for run
79  gMultiDetArray->SetParameters(run);
80  fTarget = gMultiDetArray->GetTarget();
81  fTarget->SetRandomized();
82  fIntLayer = 0;
83  fMultiLayer = (fTarget->NumberOfLayers() > 1);
84 }
85 
86 
87 
88 
91 
93 {
94  //Set projectile Z and A
95 
96  fProj.SetZandA(Z, A);
97 }
98 
99 
100 
101 
104 
106 {
107  //Set energy of projectile in MeV
108 
109  fProj.SetEnergy(e);
110  fEnergy = e;
111 }
112 
113 
114 
115 
118 
120 {
121  //Set name of detector which will detect particle
122  fDetector = gMultiDetArray->GetDetector(det);
123  if(!fDetector)
124  {
125  Error("SetDetector", "Detector %s is unknown!", det);
126  return;
127  }
128  //get list of all detectors particle must pass through to get to required detector
129  if(!fDetector->GetNode()->GetForwardTrajectories())
130  {
131  Error("SetDetector", "There is no forwards trajectory from detector %s:",det);
132  return;
133  }
134  if(fDetector->GetNode()->GetForwardTrajectories()->GetEntries()>1)
135  {
136  Error("SetDetector", "Ambiguous: there is more than one forwards trajectory from detector %s:",det);
137  fDetector->GetNode()->GetForwardTrajectories()->ls();
138  return;
139  }
140  auto traj = (KVGeoDNTrajectory*)fDetector->GetNode()->GetForwardTrajectories()->First();
141  //we store the association between detector type and index in list
142  fDetInd.Clear();
143  fAlignedDetectors.Clear();
144 
145  traj->IterateBackFrom();
146 
147  Int_t i = 0;
148  KVGeoDetectorNode* gdn = nullptr;
149  while ((gdn = traj->GetNextNode()) && gdn!=fDetector->GetNode()) {
150  KVDetector* d = gdn->GetDetector();
151  fAlignedDetectors.Add(d);
152  //check same type is not already present
153  if (fDetInd.HasParameter(d->GetType())) {
154  //detetors with same type will be called "Type_1", "Type_2" etc.
155  TString newname;
156  int j = 1;
157  newname.Form("%s_%d", d->GetType(), j++);
158  while (fDetInd.HasParameter(newname.Data()))
159  newname.Form("%s_%d", d->GetType(), j++);
160  fDetInd.SetValue(newname.Data(), i);
161  }
162  else {
163  fDetInd.SetValue(d->GetType(), i);
164  }
165  ++i;
166  }
167  //store number of aligned detectors
168  fNDets = i;
169 }
170 
171 
172 
173 
179 
181 {
182  //For multilayer targets, use this method to choose in which
183  //layer the scattering will take place.
184  //If name="", reset any previous choice so that scattering
185  //can take place in any layer
186  if (!fTarget) {
187  cout <<
188  "<KVElasticScatter::SetTargetScatteringLayer> : No target set. Set run first."
189  << endl;
190  return;
191  }
192  fIntLayer = fTarget->GetLayerIndex(name);
193  if (fIntLayer)
194  fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
195 }
196 
197 
198 
199 
203 
205 {
206  //Set binning of the GetEnergy histogram
207  // Default value is 500
208  fBinE = nbins;
209 }
210 
211 
212 
216 
218 {
219  //Perform scattering 'N' times for current values
220  //of particle Z, A and energy, target excited state, and detector.
221 
222  if (!fProj.IsDefined()) {
223  cout <<
224  "<KVElasticScatter::CalculateScattering> : Set projectile properties first"
225  << endl;
226  return;
227  }
228  if (!fEnergy) {
229  cout <<
230  "<KVElasticScatter::CalculateScattering> : Set projectile energy first"
231  << endl;
232  return;
233  }
234  if (!fDetector) {
235  cout <<
236  "<KVElasticScatter::CalculateScattering> : Set detector first" <<
237  endl;
238  return;
239  }
240  if (!fTarget) {
241  cout <<
242  "<KVElasticScatter::CalculateScattering> : No target set. Set run first."
243  << endl;
244  return;
245  }
246 
247  /* -------------------------------------------------------------------------------------------------------------------------- */
248 
249  //set up histograms
250 
251  /* -------------------------------------------------------------------------------------------------------------------------- */
252 
253  // -- histograms for debugging: depth in target and angular distribution
254  if (fDepth)
255  delete fDepth;
256  if (fTheta)
257  delete fTheta;
258  fDepth =
259  new TH1F("hDepth", "Depth (mg/cm2)", 500, 0.,
260  fTarget->GetTotalEffectiveThickness());
261  fTheta = new TH1F("hTheta", "Theta (deg.)", 500, 0., 0.);
262 
263  /* -------------------------------------------------------------------------------------------------------------------------- */
264 
265  //set up histograms for all detectors particle passes through
266  if (!fHistos) {
267  fHistos = new TObjArray(fAlignedDetectors.GetSize());
268  fHistos->SetOwner(); //will delete histos it stores
269  }
270  else {
271  fHistos->Clear(); //delete any previous histograms
272  }
273  KVDetector* d;
274  TIter n(&fAlignedDetectors);
275  while ((d = (KVDetector*) n())) {
276  fHistos->
277  Add(new
278  TH1F(Form("hEloss_%s", d->GetName()), "Eloss (MeV)", fBinE, 0., 0.));
279  }
280 
281  /* -------------------------------------------------------------------------------------------------------------------------- */
282 
283  //set up kinematics
284  if (!fKinematics)
285  fKinematics = new KV2Body;
286  fProj.SetEnergy(fEnergy);
287  fProj.SetTheta(0);
288 
289  /* -------------------------------------------------------------------------------------------------------------------------- */
290 
291  //set random interaction point for scattering
292  if (fIntLayer) {
293  fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
294  }
295  else {
296  fTarget->GetInteractionPoint(&fProj);
297  }
298 
299  /* -------------------------------------------------------------------------------------------------------------------------- */
300 
301  //get target nucleus properties from scattering layer
302  TVector3 IP = fTarget->GetInteractionPoint();
303  KVMaterial* targ_mat = fTarget->GetLayer(IP);
304  KVNucleus t;
305  t.SetZ((Int_t) targ_mat->GetZ());
306  t.SetA((Int_t) targ_mat->GetMass());
307  fKinematics->SetTarget(t);
308 
309  /* -------------------------------------------------------------------------------------------------------------------------- */
310 
311  //set excited state of target nucleus - in other words, dissipated energy for
312  //reaction due to excitation of target
313  fKinematics->SetEDiss(fExx);
314 
315  /* -------------------------------------------------------------------------------------------------------------------------- */
316 
317  Double_t xsec;
318  for (int i = 0; i < N; i++) {
319  //calculate slowing of incoming projectile
320  fTarget->SetIncoming();
321  fTarget->DetectParticle(&fProj);
322  fKinematics->SetProjectile(fProj);
323  fKinematics->SetOutgoing(fProj);
324  fKinematics->CalculateKinematics();
325  //set random direction of outgoing projectile
326 
327  double th, ph;
328  th = ph = 0.;
329  fDetector->GetRandomAngles(th, ph);
330 
331  //set energy of scattered nucleus
332  //WARNING: for inverse kinematics reactions, their are two energies for
333  //each angle below the maximum scattering angle.
334  //We only use the highest energy corresponding to the most forward CM angle.
335  Double_t e1, e2;
336  fKinematics->GetELab(3, th, 3, e1, e2);
337  fProj.SetEnergy(TMath::Max(e1, e2));
338  fProj.SetTheta(th);
339  fProj.SetPhi(ph);
340  xsec = TMath::Abs(fKinematics->GetXSecRuthLab(fProj.GetTheta()));
341  fTheta->Fill(fProj.GetTheta(), xsec);
342  //slowing of outgoing projectile in target
343  fTarget->SetOutgoing();
344  fTarget->DetectParticle(&fProj);
345  fDepth->Fill(IP.z());
346 
347  //now detect particle in detector(s)
348  int j = 0;
349  n.Reset();
350  while ((d = (KVDetector*) n())) {
351  auto cew = d->GetCentreOfEntranceWindow();
352  d->DetectParticle(&fProj, &cew);
353  //fill histograms
354  ((TH1F*)(*fHistos)[j++])->Fill(d->GetEnergy(), xsec);
355  //prepare for next round: set energy loss to zero
356  d->Clear();
357  }
358  fProj.SetEnergy(fEnergy);
359  fProj.SetTheta(0);
360  fProj.GetParameters()->Clear();
361  //set random interaction point for scattering
362  if (fIntLayer) {
363  fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
364  }
365  else {
366  fTarget->GetInteractionPoint(&fProj);
367  //if target is multilayer and the interaction layer is not fixed,
368  //the layer & hence the target nucleus may change
369  if (fMultiLayer) {
370  targ_mat = fTarget->GetLayer(fTarget->GetInteractionPoint());
371  KVNucleus t;
372  t.SetZandA((Int_t) targ_mat->GetZ(), (Int_t) targ_mat->GetMass());
373  fKinematics->SetTarget(t);
374  }
375  }
376  IP = fTarget->GetInteractionPoint();
377  }
378 }
379 
380 
381 
382 
383 
389 
391 {
392  //Energy loss in detector of given 'type' through which scattered particle passes.
393  //Warning: if there are several detectors of the same type in the list of detectors
394  //through which the particle passes, the first one (as seen by the impinging
395  //particle) will have type "type", the second "type_1", the third "type_2", etc.
396 
397  return (fDetInd.HasParameter(type) ? GetEnergy(fDetInd.GetIntValue(type)) : 0);
398 }
399 
400 
401 
402 
408 
410 {
411  //Energy loss in any detector through which scattered particle passes.
412  //The index corresponds to the order in which detectors are crossed by the
413  //particle, beginning with 0 for the first detector, and ending with
414  //(GetNDets()-1)
415  return (TH1F*)(*fHistos)[index];
416 }
417 
418 
419 //__________________________________________________________________//
420 
int Int_t
#define d(i)
#define e(i)
char Char_t
constexpr Bool_t kFALSE
double Double_t
constexpr Bool_t kTRUE
#define N
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
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 Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
char name[80]
char * Form(const char *fmt,...)
Relativistic binary kinematics calculator.
Definition: KV2Body.h:166
Calculate elastic scattering spectra in multidetector arrays.
void SetDetector(const Char_t *det)
Set name of detector which will detect particle.
TH1F * GetEnergy()
Return pointer to energy loss histogram for chosen detector (in MeV)
void CalculateScattering(Int_t N)
void SetEbinning(Int_t nbins=500)
void SetProjectile(Int_t Z, Int_t A)
Set projectile Z and A.
void SetRun(Int_t run)
Set detector parameters, target, etc. for run.
void SetTargetScatteringLayer(const Char_t *name)
void SetEnergy(Double_t E)
Set energy of projectile in MeV.
Path taken by particles through multidetector geometry.
Information on relative positions of detectors & particle trajectories.
KVDetector * GetDetector() const
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)
virtual void SetSimMode(Bool_t on=kTRUE)
Description of properties and kinematics of atomic nuclei.
Definition: KVNucleus.h:123
void SetA(Int_t a)
Definition: KVNucleus.cpp:647
void SetZ(Int_t z, Char_t mt=-1)
Definition: KVNucleus.cpp:696
void SetZandA(Int_t z, Int_t a)
Set atomic number and mass number.
Definition: KVNucleus.cpp:713
void SetRandomized(Bool_t r=kTRUE)
Definition: KVTarget.h:253
const char * Data() const
void Form(const char *fmt,...)
Double_t z() const
const Int_t n
void Error(const char *location, const char *fmt,...)
void Warning(const char *location, const char *fmt,...)
Add
Double_t Abs(Double_t d)
Double_t Max(Double_t a, Double_t b)
ClassImp(TPyArg)