KaliVeda
Toolkit for HIC analysis
KVDetector.cpp
1 #include "KVDetector.h"
2 #include "TROOT.h"
3 #include "KVGroup.h"
4 #include "KVCalibrator.h"
5 #include "TPluginManager.h"
6 #include "TClass.h"
7 /*** geometry ***/
8 #include "TGeoVolume.h"
9 #include "TGeoManager.h"
10 #include "TGeoMatrix.h"
11 #include "TGeoBBox.h"
12 #include "TGeoArb8.h"
13 #include "TGeoTube.h"
14 #include "KVCalibratedSignal.h"
15 #include "KVDetectorSignalExpression.h"
16 #include "KVZDependentCalibratedSignal.h"
17 #include "KVMaterialStack.h"
18 #include <TGeoPhysicalNode.h>
19 #include <TGraph.h>
20 
21 using namespace std;
22 
23 ClassImp(KVDetector)
24 
25 Int_t KVDetector::fDetCounter = 0;
26 
27 
30 
31 void KVDetector::init()
32 {
33  //default initialisations
34  fCalibrators = nullptr;
35  fParticles = nullptr;
36  fActiveLayer = nullptr;
37  fIdentP = fUnidentP = 0;
38  fELossF = fEResF = fRangeF = nullptr;
39  fEResforEinc = -1.;
40  fSimMode = kFALSE;
41  fPresent = kTRUE;
42  fDetecting = kTRUE;
43  fParentStrucList.SetCleanup();
44  fSingleLayer = kTRUE;
45  fNode.SetDetector(this);
46  // detector owns any signals which are defined for it
47  fDetSignals.SetOwner();
48  // adding a new signal with the same name as an existing one
49  // will delete the existing signal and replace it
50  fDetSignals.ReplaceObjects();
51 }
52 
53 
54 
57 
58 KVDetector::KVDetector()
59 {
60 //default ctor
61  init();
62  fDetCounter++;
63  SetName(Form("Det_%d", fDetCounter));
64 }
65 
66 
67 
70 
71 KVDetector::KVDetector(const Char_t* type,
72  const Float_t thick): KVMaterial()
73 {
74  // Create a new detector of a given material and thickness in centimetres (default value = 0.0)
75 
76  init();
77  SetType("DET");
78  fDetCounter++;
79  SetName(Form("Det_%d", fDetCounter));
80  AddAbsorber(new KVMaterial(type, thick));
81 }
82 
83 
84 
89 
90 KVDetector::KVDetector(const Char_t* gas, const Double_t thick, const Double_t pressure, const Double_t temperature)
91 {
92  // Create gaseous dteector with given type, linear thickness in cm, pressure in Torr, and temperature in degrees C (default value 19°C).
93  //
94  // \note this just defines some gas, the 'detector' has no windows!
95 
96  init();
97  SetType("DET");
98  fDetCounter++;
99  SetName(Form("Det_%d", fDetCounter));
100  AddAbsorber(new KVMaterial(gas, thick, pressure, temperature));
101 }
102 
103 
104 
107 
108 KVDetector::KVDetector(const KVDetector& obj) : KVMaterial(), KVPosition()
109 {
110 //copy ctor
111  init();
112  fDetCounter++;
113  SetName(Form("Det_%d", fDetCounter));
114 #if ROOT_VERSION_CODE >= ROOT_VERSION(3,4,0)
115  obj.Copy(*this);
116 #else
117  ((KVDetector&) obj).Copy(*this);
118 #endif
119 }
120 
121 
122 #if ROOT_VERSION_CODE >= ROOT_VERSION(3,4,0)
123 
128 
129 void KVDetector::Copy(TObject& obj) const
130 #else
131 void KVDetector::Copy(TObject& obj)
132 #endif
133 {
134  //copy 'this' to 'obj'
135  //The structure of the detector is copied, with new cloned objects for
136  //each absorber layer. The active layer is set in the new detector.
137 
138  TIter next(&fAbsorbers);
139  KVMaterial* mat;
140  while ((mat = (KVMaterial*) next())) {
141  ((KVDetector&) obj).AddAbsorber((KVMaterial*) mat->Clone());
142  }
143  //set active layer
144  Int_t in_actif = fAbsorbers.IndexOf(fActiveLayer);
145  ((KVDetector&) obj).SetActiveLayer(((KVDetector&)obj).GetAbsorber(in_actif));
146 }
147 
148 
149 
150 
152 
153 KVDetector::~KVDetector()
154 {
155  SafeDelete(fCalibrators);
156  SafeDelete(fParticles);
157  SafeDelete(fELossF);
158  SafeDelete(fEResF);
159  SafeDelete(fRangeF);
160 }
161 
162 
163 
168 
169 void KVDetector::SetMaterial(const Char_t* type)
170 {
171  // Set material of active layer.
172  // If no absorbers have been added to the detector, create and add
173  // one (active layer by default)
174 
175  if (!GetActiveLayer())
176  AddAbsorber(new KVMaterial(type));
177  else
178  GetActiveLayer()->SetMaterial(type);
179 }
180 
181 
182 
194 
195 void KVDetector::DetectParticle(KVNucleus* kvp, TVector3* norm)
196 {
197  //Calculate the energy loss of a charged particle traversing the detector,
198  //the particle is slowed down, it is added to the list of all particles hitting the
199  //detector. The apparent energy loss of the particle in the active layer of the
200  //detector is set.
201  //Do nothing if particle has zero (or -ve) energy.
202  //
203  //If the optional argument 'norm' is given, it is supposed to be a vector
204  //normal to the detector, oriented from the origin towards the detector.
205  //In this case the effective thicknesses of the detector's absorbers 'seen' by the particle
206  //depending on its direction of motion is used for the calculation.
207 
208  if (kvp->GetKE() <= 0.)
209  return;
210 
211  AddHit(kvp); //add nucleus to list of particles hitting detector in the event
212  //set flag to say that particle has been slowed down
213  kvp->SetIsDetected();
214  //If this is the first absorber that the particle crosses, we set a "reminder" of its
215  //initial energy
216  if (!kvp->GetPInitial())
217  kvp->SetE0();
218 
219  std::vector<Double_t> thickness;
220  if (norm) {
221  // modify thicknesses of all layers according to orientation,
222  // and store original thicknesses in array
223  TVector3 p = kvp->GetMomentum();
224  thickness.reserve(fAbsorbers.GetEntries());
225  KVMaterial* abs;
226  int i = 0;
227  TIter next(&fAbsorbers);
228  while ((abs = (KVMaterial*) next())) {
229  thickness[i++] = abs->GetThickness();
230  Double_t T = abs->GetEffectiveThickness((*norm), p);
231  abs->SetThickness(T);
232  }
233  }
234  Double_t eloss = GetTotalDeltaE(kvp->GetZ(), kvp->GetA(), kvp->GetEnergy());
235  Double_t dE = GetDeltaE(kvp->GetZ(), kvp->GetA(), kvp->GetEnergy());
236  if (norm) {
237  // reset thicknesses of absorbers
238  KVMaterial* abs;
239  int i = 0;
240  TIter next(&fAbsorbers);
241  while ((abs = (KVMaterial*) next())) {
242  abs->SetThickness(thickness[i++]);
243  }
244  }
245  Double_t epart = kvp->GetEnergy() - eloss;
246  if (epart < 1e-3) {
247  //printf("%s, pb d arrondi on met l energie de la particule a 0\n",GetName());
248  epart = 0.0;
249  }
250  kvp->SetEnergy(epart);
251  Double_t eloss_old = GetEnergyLoss();
252  SetEnergyLoss(eloss_old + dE);
253 }
254 
255 
256 
257 
267 
268 Double_t KVDetector::GetELostByParticle(KVNucleus* kvp, TVector3* norm)
269 {
270  //Calculate the total energy loss of a charged particle traversing the detector.
271  //This does not affect the "stored" energy loss value of the detector,
272  //nor its ACQData, nor the energy of the particle.
273  //
274  //If the optional argument 'norm' is given, it is supposed to be a vector
275  //normal to the detector, oriented from the origin towards the detector.
276  //In this case the effective thicknesses of the detector's absorbers 'seen' by the particle
277  //depending on its direction of motion is used for the calculation.
278 
279  std::vector<Double_t> thickness;
280  if (norm) {
281  // modify thicknesses of all layers according to orientation,
282  // and store original thicknesses in array
283  TVector3 p = kvp->GetMomentum();
284  thickness.reserve(fAbsorbers.GetEntries());
285  KVMaterial* abs;
286  int i = 0;
287  TIter next(&fAbsorbers);
288  while ((abs = (KVMaterial*) next())) {
289  thickness[i++] = abs->GetThickness();
290  Double_t T = abs->GetEffectiveThickness((*norm), p);
291  abs->SetThickness(T);
292  }
293  }
294  Double_t eloss = GetTotalDeltaE(kvp->GetZ(), kvp->GetA(), kvp->GetEnergy());
295  if (norm) {
296  // reset thicknesses of absorbers
297  KVMaterial* abs;
298  int i = 0;
299  TIter next(&fAbsorbers);
300  while ((abs = (KVMaterial*) next())) {
301  abs->SetThickness(thickness[i++]);
302  }
303  }
304  return eloss;
305 }
306 
307 
308 
309 
319 
320 Double_t KVDetector::GetParticleEIncFromERes(KVNucleus* kvp, TVector3* norm)
321 {
322  //Calculate the energy of particle 'kvn' before its passage through the detector,
323  //based on the current kinetic energy, Z & A of nucleus 'kvn', supposed to be
324  //after passing through the detector.
325  //
326  //If the optional argument 'norm' is given, it is supposed to be a vector
327  //normal to the detector, oriented from the origin towards the detector.
328  //In this case the effective thicknesses of the detector's absorbers 'seen' by the particle
329  //depending on its direction of motion is used for the calculation.
330 
331  Double_t Einc = 0.0;
332  //make 'clone' of particle
333  KVNucleus clone_part(kvp->GetZ(), kvp->GetA());
334  clone_part.SetMomentum(kvp->GetMomentum());
335  //composite detector - calculate losses in all layers
336  KVMaterial* abs;
337  TIter next(&fAbsorbers, kIterBackward); //work through list backwards
338  while ((abs = (KVMaterial*) next())) {
339 
340  //calculate energy of particle before current absorber
341  Einc = abs->GetParticleEIncFromERes(&clone_part, norm);
342 
343  //set new energy of particle
344  clone_part.SetKE(Einc);
345  }
346  return Einc;
347 }
348 
349 
350 
351 
355 
356 void KVDetector::Print(Option_t* opt) const
357 {
358  //Print info on this detector
359  //if option="data" the energy loss and raw data are displayed
360 
361  if (!strcmp(opt, "data")) {
362  cout << GetName() << " -- ";
363  // print raw data signals
364  KVDetectorSignal* sig;
365  TIter it(&GetListOfDetectorSignals());
366  while ((sig = (KVDetectorSignal*)it())) {
367  if (sig->IsRaw()) {
368  std::cout << sig->GetName() << "=" << sig->GetValue() << " | ";
369  }
370  }
371  // print calibrated data
372  it.Reset();
373  while ((sig = (KVDetectorSignal*)it())) {
374  if (!sig->IsRaw() && !sig->GetValueNeedsExtraParameters()) {
375  std::cout << sig->GetName() << "=" << sig->GetValue() << " | ";
376  }
377  }
378  cout << " ";
379  if (BelongsToUnidentifiedParticle())
380  cout << "(Belongs to an unidentified particle)";
381  cout << endl;
382  }
383  else if (!strcmp(opt, "all")) {
384  //Give full details of detector
385  //
386  TString option(" ");
387  cout << option << ClassName() << " : " << ((KVDetector*) this)->
388  GetName() << endl;
389  //composite detector - print all layers
390  KVMaterial* abs;
391  TIter next(&fAbsorbers);
392  while ((abs = (KVMaterial*) next())) {
393  if (GetActiveLayer() == abs)
394  cout << " ### ACTIVE LAYER ### " << endl;
395  cout << option << option;
396  abs->Print();
397  if (GetActiveLayer() == abs)
398  cout << " #################### " << endl;
399  }
400  if (fParticles) {
401  cout << option << " --- Detected particles:" << endl;
402  fParticles->Print();
403  }
404  }
405  else {
406  //just print name
407  cout << ClassName() << ": " << ((KVDetector*) this)->
408  GetName() << endl;
409  }
410 }
411 
412 
413 
438 
439 Bool_t KVDetector::AddCalibrator(KVCalibrator* cal, const KVNameValueList& opts)
440 {
441  // Associate a calibration with this detector.
442  //
443  // This will add a new signal to the list of the detector's signals.
444  //
445  // Also sets calibrator's name to `[detname]_[caltype]` where `caltype` is the type of the KVCalibration object.
446  //
447  // \param[in] cal pointer to KVCalibrator object (must be on heap, i.e. created with new: detector handles deletion)
448  // \param[in] opts
449  // \parblock
450  // can be used to pass any extra parameters/options needed by the calibrator.
451  // For example, if it contains a parameter `ZRange`:
452  //
453  //~~~~~~~~~~~~~~~
454  // ZRange=1-10
455  //~~~~~~~~~~~~~~~
456  //
457  // then the calibrator will be handled by a KVZDependentCalibratedSignal (handles several
458  // calibrators which provide the same output signal, each one is used for a specific range
459  // of atomic numbers)
460  // \endparblock
461  //
462  // \returns kFALSE in case of problems (non-existent input signal for calibrator, output signal not defined for calibrator),
463  // otherwise kTRUE
464 
465  if (!cal) return kFALSE;
466  // look for input signal
467  KVDetectorSignal* in = GetDetectorSignal(cal->GetInputSignalType());
468  // input signal must exist
469  if (!in) {
470  Error("AddCalibrator", "Detector:%s has no signal %s, required as input by calibrator %s to provide output %s",
471  GetName(), cal->GetInputSignalType().Data(), cal->GetType(), cal->GetOutputSignalType().Data());
472  return kFALSE;
473  }
474  if (!fCalibrators)
475  fCalibrators = new KVList();
476 
477  cal->SetDetector(this);
478  fCalibrators->Add(cal);
479  cal->SetName(Form("%s_%s", GetName(), cal->GetType()));
480 
481  if (cal->GetOutputSignalType() == "") {
482  Warning("AddCalibrator", "%s : output signal not defined for calibrator %s. No output signal created.",
483  GetName(), cal->GetType());
484  }
485  else {
486  KVDetectorSignal* new_cal_sig(nullptr);
487  if (opts.HasParameter("ZRange")) {
488  // If 'ZRange' parameter is given we need to find the KVZDependentCalibratedSignal
489  // and add a new signal to it.
490  KVDetectorSignal* sig = GetDetectorSignal(cal->GetOutputSignalType());
491  if (!sig) {
492  new_cal_sig = sig = new KVZDependentCalibratedSignal(in, cal->GetOutputSignalType());
493  }
494  dynamic_cast<KVZDependentCalibratedSignal*>(sig)->AddSignal(new KVCalibratedSignal(in, cal), opts.GetStringValue("ZRange"));
495  }
496  else {
497  new_cal_sig = new KVCalibratedSignal(in, cal);
498  }
499  if (new_cal_sig) fDetSignals.Add(new_cal_sig);
500  }
501  return kTRUE;
502 }
503 
504 
505 
516 
517 Bool_t KVDetector::ReplaceCalibrator(const Char_t* type, KVCalibrator* cal, const KVNameValueList& opts)
518 {
519  // Replace calibrator of given type with the given calibrator object
520  // The calibrator object should not be shared with any other detectors: it now belongs
521  // to this detector, which will delete it when necessary.
522  // If an exising calibrator with the same type is already defined, it will be
523  // deleted and removed from the detector's calibrator list
524  //
525  // Returns kFALSE in case of problems.
526  //
527  // The (optional) KVNameValueList argument can be used to pass any extra parameters/options.
528 
529  KVCalibrator* old_cal = GetCalibrator(type);
530  if (old_cal) {
531  fCalibrators->Remove(old_cal);
532  remove_signal_for_calibrator(old_cal);
533  delete old_cal;
534  }
535  return AddCalibrator(cal, opts);
536 }
537 
538 
539 
544 
545 Bool_t KVDetector::IsCalibrated(const KVNameValueList& params) const
546 {
547  // A detector is considered to be calibrated if it has
548  // a signal "Energy" available and if depending on the supplied parameters
549  // this signal can be calculated
550 
551  KVCalibratedSignal* e_sig = dynamic_cast<KVCalibratedSignal*>(GetDetectorSignal("Energy"));
552  return (e_sig && e_sig->IsAvailableFor(params) && IsOK());
553 }
554 
555 
556 
557 
564 
565 void KVDetector::Clear(Option_t* opt)
566 {
567  //Set energy loss(es) etc. to zero
568  //
569  //If opt="N":
570  // + we do not reset acquisition parameters/raw detector signals
571  // + in SimMode we do not reset energy losses (this is called before reconstruction happens)
572 
573  bool option_N = (strncmp(opt, "N", 1) == 0);
574  bool sim_mode_option_N = option_N && IsSimMode();
575 
576  if (!sim_mode_option_N) KVMaterial::Clear(opt);
577  SetAnalysed(kFALSE);
578  fIdentP = fUnidentP = 0;
579  ResetBit(kIdentifiedParticle);
580  ResetBit(kUnidentifiedParticle);
581  if (!option_N) {
582  TIter it(&GetListOfDetectorSignals());
583  KVDetectorSignal* ds;
584  while ((ds = (KVDetectorSignal*)it())) {
585  if (ds->IsRaw() && !ds->IsExpression()) ds->Reset();
586  }
587  }
588  ClearHits();
589  if (!sim_mode_option_N) {
590  //reset all layers in detector
591  KVMaterial* mat;
592  TIter next(&fAbsorbers);
593  while ((mat = (KVMaterial*) next())) {
594  mat->Clear();
595  }
596  }
597  fEResforEinc = -1.;
598 }
599 
600 
601 
606 
607 void KVDetector::AddAbsorber(KVMaterial* mat)
608 {
609  // Add a layer of absorber material to the detector
610  // By default, the first layer added is set as the "Active" layer.
611  // Call SetActiveLayer to change this.
612  fAbsorbers.Add(mat);
613  if (!TestBit(kActiveSet))
614  SetActiveLayer(mat);
615  if (fAbsorbers.GetSize() > 1) fSingleLayer = kFALSE;
616 }
617 
618 
619 
622 
623 KVMaterial* KVDetector::GetAbsorber(Int_t i) const
624 {
625  //Returns pointer to the i-th absorber in the detector (i=0 first absorber, i=1 second, etc.)
626 
627  if (!fAbsorbers.GetEntries()) {
628  Error("GetAbsorber", "No absorbers defined for detector");
629  return nullptr;
630  }
631  if (i < 0 || i >= fAbsorbers.GetEntries()) {
632  Error("GetAbsorber", "i=%d but valid values are 0-%d [number of absorbers in list]", i,
633  fAbsorbers.GetEntries());
634  return nullptr;
635  }
636  return (KVMaterial*) fAbsorbers.At(i);
637 }
638 
639 
640 
645 
646 void KVDetector::RemoveAllAbsorbers()
647 {
648  // Completely reset the KVDetector as if it had just been created by a call to the default constructor:
649  // + removes all absorber layers
650  // + removes all calibrators/detector signals
651 
652  Clear();
653  RemoveCalibrators();
654  SafeDelete(fCalibrators);
655  SafeDelete(fParticles);
656  fAbsorbers.Clear();
657  SafeDelete(fELossF);
658  SafeDelete(fEResF);
659  SafeDelete(fRangeF);
660  fCalibrators = nullptr;
661  fParticles = nullptr;
662  fActiveLayer = nullptr;
663  ResetBit(kActiveSet);
664  fIdentP = fUnidentP = 0;
665  fELossF = fEResF = fRangeF = nullptr;
666  fEResforEinc = -1.;
667  fSimMode = kFALSE;
668  fPresent = kTRUE;
669  fDetecting = kTRUE;
670  fSingleLayer = kTRUE;
671 }
672 
673 
674 
678 
679 void KVDetector::remove_signal_for_calibrator(KVCalibrator* K)
680 {
681  // Used when a calibrator object is removed or replaced
682  // We remove and delete the corresponding output signal from the list of detector signals
683 
684  if (K->GetOutputSignalType() != "") {
685  KVDetectorSignal* ds = GetDetectorSignal(K->GetOutputSignalType());
686  if (ds) {
687  fDetSignals.Remove(ds);
688  delete ds;
689  }
690  }
691 }
692 
693 
694 
700 
701 void KVDetector::RemoveCalibrators()
702 {
703  // Removes all calibrations associated to this detector: in other words, we delete all
704  // the KVCalibrator objects in list fCalibrators.
705  //
706  // We also destroy all signals provided by these calibrators
707 
708  if (fCalibrators) {
709  KVCalibrator* K;
710  TIter it(fCalibrators);
711  while ((K = (KVCalibrator*)it())) {
712  remove_signal_for_calibrator(K);
713  }
714  fCalibrators->Delete();
715  }
716 }
717 
718 
719 
720 
750 
751 Double_t KVDetector::GetCorrectedEnergy(KVNucleus* nuc, Double_t e, Bool_t transmission)
752 {
753  // Returns the total energy loss in the detector for a given nucleus
754  // including inactive absorber layers.
755  // e = energy loss in active layer (if not given, we use current value)
756  // transmission = kTRUE (default): the particle is assumed to emerge with
757  // a non-zero residual energy Eres after the detector.
758  // = kFALSE: the particle is assumed to stop in the detector.
759  //
760  // WARNING: if transmission=kTRUE, and if the residual energy after the detector
761  // is known (i.e. measured in a detector placed after this one), you should
762  // first call
763  // SetEResAfterDetector(Eres);
764  // before calling this method. Otherwise, especially for heavy ions, the
765  // correction may be false for particles which are just above the punch-through energy.
766  //
767  // WARNING 2: if measured energy loss in detector active layer is greater than
768  // maximum possible theoretical value for given nucleus' Z & A, this may be because
769  // the A was not measured but calculated from Z and hence could be false, or perhaps
770  // there was an (undetected) pile-up of two or more particles in the detector.
771  // In this case we return the corrected energy
772  // corresponding to the maximum theoretical energy loss in the active layer
773  // and we add the following parameters to the particle (in nuc->GetParameters()):
774  //
775  // GetCorrectedEnergy.Warning = 1
776  // GetCorrectedEnergy.Detector = [name]
777  // GetCorrectedEnergy.MeasuredDE = [value]
778  // GetCorrectedEnergy.MaxDE = [value]
779  // GetCorrectedEnergy.Transmission = 0 or 1
780  // GetCorrectedEnergy.ERES = [value]
781 
782  Int_t z = nuc->GetZ();
783  Int_t a = nuc->GetA();
784 
785  if (e < 0.) e = GetEnergy();
786  if (e <= 0) {
787  SetEResAfterDetector(-1.);
788  return 0;
789  }
790 
791  // check if calling this function for a previous detector in the particle's trajectory
792  // led to using an effective incident angle for the particle: if so, use it here
793  if (nuc->GetParameters()->HasDoubleParameter("GetCorrectedEnergy.IncidentAngle")) {
794  KVMaterialStack stack(this);
795  stack.SetIncidenceAngle(nuc->GetParameters()->GetDoubleValue("GetCorrectedEnergy.IncidentAngle"));
796  // check that apparent energy loss in detector is compatible with a & z
797  Double_t maxDE = stack.GetMaxDeltaE(z, a);
798  if (e > maxDE) {
799  auto inc_angle = stack.GetMinimumIncidentAngleForDEMax(z, a, e);
800  stack.SetIncidenceAngle(inc_angle);
801  nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.Warning", GetName()), 1);
802  nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.MeasuredDE", GetName()), e);
803  nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.MaxDE", GetName()), maxDE);
804  nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.Transmission", GetName()), (Int_t)transmission);
805  nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.ERES", GetName()), GetEResAfterDetector());
806  nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.IncidentAngle", GetName()), inc_angle);
807  nuc->GetParameters()->SetValue("GetCorrectedEnergy.IncidentAngle", inc_angle);
808  }
809  return get_corrected_energy(&stack, nuc, e, transmission);
810  }
811  // check that apparent energy loss in detector is compatible with a & z
812  Double_t maxDE = GetMaxDeltaE(z, a);
813  if (e > maxDE) {
814  KVMaterialStack stack(this);
815  auto inc_angle = stack.GetMinimumIncidentAngleForDEMax(z, a, e);
816  stack.SetIncidenceAngle(inc_angle);
817  nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.Warning", GetName()), 1);
818  nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.MeasuredDE", GetName()), e);
819  nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.MaxDE", GetName()), maxDE);
820  nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.Transmission", GetName()), (Int_t)transmission);
821  nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.ERES", GetName()), GetEResAfterDetector());
822  nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.IncidentAngle", GetName()), inc_angle);
823  nuc->GetParameters()->SetValue("GetCorrectedEnergy.IncidentAngle", inc_angle);
824  return get_corrected_energy(&stack, nuc, e, transmission);
825  }
826  return get_corrected_energy(this, nuc, e, transmission);
827 }
828 
829 
830 
831 
849 
850 Int_t KVDetector::FindZmin(Double_t ELOSS, Char_t mass_formula)
851 {
852  //For particles which stop in the first stage of an identification telescope,
853  //we can at least estimate a minimum Z value based on the energy lost in this
854  //detector.
855  //
856  //This is based on the KVMaterial::GetMaxDeltaE method, giving the maximum
857  //energy loss in the active layer of the detector for a given nucleus (A,Z).
858  //
859  //The "Zmin" is the Z of the particle which gives a maximum energy loss just greater
860  //than that measured in the detector. Particles with Z<Zmin could not lose as much
861  //energy and so are excluded.
862  //
863  //If ELOSS is not given, we use the current value of GetEnergy()
864  //Use 'mass_formula' to change the formula used to calculate the A of the nucleus
865  //from its Z. Default is valley of stability value. (see KVNucleus::GetAFromZ).
866  //
867  //If the value of ELOSS or GetEnergy() is <=0 we return Zmin=0
868 
869  ELOSS = (ELOSS < 0. ? GetEnergy() : ELOSS);
870  if (ELOSS <= 0) return 0;
871 
872  UInt_t z = 40;
873  UInt_t zmin, zmax;
874  zmin = 1;
875  zmax = 92;
876  Double_t difference;
877  UInt_t last_positive_difference_z = 1;
878  KVNucleus particle;
879  if (mass_formula > -1)
880  particle.SetMassFormula((UChar_t)mass_formula);
881 
882  difference = 0.;
883 
884  while (zmax > zmin + 1) {
885 
886  particle.SetZ(z);
887 
888  difference = GetMaxDeltaE(z, particle.GetA()) - ELOSS;
889  //if difference < 0 the z is too small
890  if (difference < 0.0) {
891 
892  zmin = z;
893  z += (UInt_t)((zmax - z) / 2 + 0.5);
894 
895  }
896  else {
897 
898  zmax = z;
899  last_positive_difference_z = z;
900  z -= (UInt_t)((z - zmin) / 2 + 0.5);
901 
902  }
903  }
904  if (difference < 0) {
905  //if the last z we tried actually made the difference become negative again
906  //(i.e. z too small) we return the last z which gave a +ve difference
907  z = last_positive_difference_z;
908  }
909  return z;
910 }
911 
912 
913 
914 
923 
924 Double_t KVDetector::ELossActive(Double_t* x, Double_t* par)
925 {
926  // Calculates energy loss (in MeV) in active layer of detector, taking into account preceding layers
927  //
928  // Arguments are:
929  // x[0] is incident energy in MeV
930  // Parameters are:
931  // par[0] Z of ion
932  // par[1] A of ion
933 
934  Double_t e = x[0];
935  TIter next(&fAbsorbers);
936  KVMaterial* mat;
937  mat = (KVMaterial*)next();
938  while (fActiveLayer != mat) {
939  // calculate energy losses in absorbers before active layer
940  e = mat->GetERes(par[0], par[1], e); //residual energy after layer
941  if (e <= 0.)
942  return 0.; // return 0 if particle stops in layers before active layer
943  mat = (KVMaterial*)next();
944  }
945  //calculate energy loss in active layer
946  return mat->GetDeltaE(par[0], par[1], e);
947 }
948 
949 
950 
951 
961 
962 Double_t KVDetector::RangeDet(Double_t* x, Double_t* par)
963 {
964  // Calculates range (in centimetres) of ions in detector as a function of incident energy (in MeV),
965  // taking into account all layers of the detector.
966  //
967  // Arguments are:
968  // x[0] = incident energy in MeV
969  // Parameters are:
970  // par[0] = Z of ion
971  // par[1] = A of ion
972 
973  Double_t Einc = x[0];
974  Int_t Z = (Int_t)par[0];
975  Int_t A = (Int_t)par[1];
976  Double_t range = 0.;
977  TIter next(&fAbsorbers);
978  KVMaterial* mat = (KVMaterial*)next();
979  if (!mat) return 0.;
980  do {
981  // range in this layer
982  Double_t this_range = mat->GetLinearRange(Z, A, Einc);
983  KVMaterial* next_mat = (KVMaterial*)next();
984  if (this_range > mat->GetThickness()) {
985  // particle traverses layer.
986  if (next_mat)
987  range += mat->GetThickness();
988  else // if this is last layer, the range continues to increase beyond size of detector
989  range += this_range;
990  // calculate incident energy for next layer (i.e. residual energy after this layer)
991  Einc = mat->GetERes(Z, A, Einc);
992  }
993  else {
994  // particle stops in layer
995  range += this_range;
996  return range;
997  }
998  mat = next_mat;
999  }
1000  while (mat);
1001  // particle traverses detector
1002  return range;
1003 }
1004 
1005 
1006 
1007 
1017 
1018 Double_t KVDetector::EResDet(Double_t* x, Double_t* par)
1019 {
1020  // Calculates residual energy (in MeV) of particle after traversing all layers of detector.
1021  // Returned value is -1000 if particle stops in one of the layers of the detector.
1022  //
1023  // Arguments are:
1024  // x[0] is incident energy in MeV
1025  // Parameters are:
1026  // par[0] Z of ion
1027  // par[1] A of ion
1028 
1029  Double_t e = x[0];
1030  TIter next(&fAbsorbers);
1031  KVMaterial* mat;
1032  while ((mat = (KVMaterial*)next())) {
1033  Double_t eres = mat->GetERes(par[0], par[1], e); //residual energy after layer
1034  if (eres <= 0.)
1035  return -1000.; // return -1000 if particle stops in layers before active layer
1036  e = eres;
1037  }
1038  return e;
1039 }
1040 
1041 
1042 
1043 
1058 
1059 KVDetector* KVDetector::MakeDetector(const Char_t* name, Float_t thickness)
1060 {
1061  //Static function which will create an instance of the KVDetector-derived class
1062  //corresponding to 'name'
1063  //These are defined as 'Plugin' objects in the file $KVROOT/KVFiles/.kvrootrc :
1064  // [name_of_dataset].detector_type
1065  // detector_type
1066  // To use the dataset-dependent plugin, call this method with
1067  // name = "[name_of_dataset].detector_type"
1068  // If not, the default plugin will be used
1069  //first we check if there is a special plugin for the DataSet
1070  //if not we take the default one
1071  //
1072  //'thickness' is passed as argument to the constructor for the detector plugin
1073 
1074  //check and load plugin library
1075  TPluginHandler* ph = nullptr;
1076  KVString nom(name);
1077  if (!nom.Contains(".") && !(ph = LoadPlugin("KVDetector", name))) return nullptr;
1078  if (nom.Contains(".")) {
1079  // name format like [dataset].[det_type]
1080  // in case dataset name contains "." we parse string to find detector type: assumed after the last "."
1081  nom.RBegin(".");
1082  KVString det_type = nom.RNext();
1083  if (!(ph = LoadPlugin("KVDetector", name))) {
1084  if (!(ph = LoadPlugin("KVDetector", det_type))) {
1085  return nullptr;
1086  }
1087  }
1088  }
1089 
1090  //execute constructor/macro for detector
1091  return (KVDetector*) ph->ExecPlugin(1, thickness);
1092 }
1093 
1094 
1095 
1098 
1099 Double_t KVDetector::GetEntranceWindowSurfaceArea()
1100 {
1101  // Return surface area of first layer of detector in cm2.
1102 
1103  if (!GetEntranceWindow().GetShape()->InheritsFrom("TGeoArb8")) {
1104  // simple shape area
1105  return GetEntranceWindow().GetShape()->GetFacetArea(1);
1106  }
1107  // Monte Carlo calculation for TGeoArb8 shapes
1108  return GetEntranceWindow().GetSurfaceArea();
1109 }
1110 
1111 
1112 
1117 
1118 TF1* KVDetector::GetEResFunction(Int_t Z, Int_t A)
1119 {
1120  // Return pointer toTF1 giving residual energy after detector as function of incident energy,
1121  // for a given nucleus (Z,A).
1122  // The TF1::fNpx parameter is taken from environment variable KVDetector.ResidualEnergy.Npx
1123 
1124  if (!fEResF) {
1125  fEResF = new TF1(Form("KVDetector:%s:ERes", GetName()), this, &KVDetector::EResDet,
1126  0., 1.e+04, 2, "KVDetector", "EResDet");
1127  fEResF->SetNpx(gEnv->GetValue("KVDetector.ResidualEnergy.Npx", 20));
1128  }
1129  fEResF->SetParameters((Double_t)Z, (Double_t)A);
1130  fEResF->SetRange(0., GetSmallestEmaxValid(Z, A));
1131  fEResF->SetTitle(Form("Residual energy [MeV] after detector %s for Z=%d A=%d", GetName(), Z, A));
1132 
1133  return fEResF;
1134 }
1135 
1136 
1137 
1142 
1143 TF1* KVDetector::GetRangeFunction(Int_t Z, Int_t A)
1144 {
1145  // Return pointer toTF1 giving range (in centimetres) in detector as function of incident energy,
1146  // for a given nucleus (Z,A).
1147  // The TF1::fNpx parameter is taken from environment variable KVDetector.Range.Npx
1148 
1149  if (!fRangeF) {
1150  fRangeF = new TF1(Form("KVDetector:%s:Range", GetName()), this, &KVDetector::RangeDet,
1151  0., 1.e+04, 2, "KVDetector", "RangeDet");
1152  fRangeF->SetNpx(gEnv->GetValue("KVDetector.Range.Npx", 20));
1153  }
1154  fRangeF->SetParameters((Double_t)Z, (Double_t)A);
1155  fRangeF->SetRange(0., GetSmallestEmaxValid(Z, A));
1156  fRangeF->SetTitle(Form("Range [cm] in detector %s for Z=%d A=%d", GetName(), Z, A));
1157 
1158  return fRangeF;
1159 }
1160 
1161 
1162 
1167 
1168 TF1* KVDetector::GetELossFunction(Int_t Z, Int_t A)
1169 {
1170  // Return pointer to TF1 giving energy loss in active layer of detector as function of incident energy,
1171  // for a given nucleus (Z,A).
1172  // The TF1::fNpx parameter is taken from environment variable KVDetector.EnergyLoss.Npx
1173 
1174  if (!fELossF) {
1175  fELossF = new TF1(Form("KVDetector:%s:ELossActive", GetName()), this, &KVDetector::ELossActive,
1176  0., 1.e+04, 2, "KVDetector", "ELossActive");
1177  fELossF->SetNpx(gEnv->GetValue("KVDetector.EnergyLoss.Npx", 20));
1178  }
1179  fELossF->SetParameters((Double_t)Z, (Double_t)A);
1180  fELossF->SetRange(0., GetSmallestEmaxValid(Z, A));
1181  fELossF->SetTitle(Form("Energy loss [MeV] in detector %s for Z=%d A=%d", GetName(), Z, A));
1182  return fELossF;
1183 }
1184 
1185 
1186 
1191 
1192 Double_t KVDetector::GetEIncOfMaxDeltaE(Int_t Z, Int_t A)
1193 {
1194  // Overrides KVMaterial::GetEIncOfMaxDeltaE
1195  // Returns incident energy corresponding to maximum energy loss in the
1196  // active layer of the detector, for a given nucleus.
1197 
1198  return GetELossFunction(Z, A)->GetMaximumX();
1199 }
1200 
1201 
1202 
1207 
1208 Double_t KVDetector::GetMaxDeltaE(Int_t Z, Int_t A)
1209 {
1210  // Overrides KVMaterial::GetMaxDeltaE
1211  // Returns maximum energy loss in the
1212  // active layer of the detector, for a given nucleus.
1213 
1214  return GetELossFunction(Z, A)->GetMaximum();
1215 }
1216 
1217 
1218 
1223 
1224 Double_t KVDetector::GetDeltaE(Int_t Z, Int_t A, Double_t Einc, Double_t)
1225 {
1226  // Overrides KVMaterial::GetDeltaE
1227  // Returns energy loss of given nucleus in the active layer of the detector.
1228 
1229  // optimization for single-layer detectors
1230  if (fSingleLayer) {
1231  return fActiveLayer->GetDeltaE(Z, A, Einc);
1232  }
1233  return GetELossFunction(Z, A)->Eval(Einc);
1234 }
1235 
1236 
1237 
1241 
1242 Double_t KVDetector::GetTotalDeltaE(Int_t Z, Int_t A, Double_t Einc)
1243 {
1244  // Returns calculated total energy loss of ion in ALL layers of the detector.
1245  // This is just (Einc - GetERes(Z,A,Einc))
1246 
1247  return Einc - GetERes(Z, A, Einc);
1248 }
1249 
1250 
1251 
1256 
1257 Double_t KVDetector::GetERes(Int_t Z, Int_t A, Double_t Einc, Double_t)
1258 {
1259  // Overrides KVMaterial::GetERes
1260  // Returns residual energy of given nucleus after the detector.
1261  // Returns 0 if Einc<=0
1262 
1263  if (Einc <= 0.) return 0.;
1264  Double_t eres = GetEResFunction(Z, A)->Eval(Einc);
1265  // Eres function returns -1000 when particle stops in detector,
1266  // in order for function inversion (GetEIncFromEres) to work
1267  if (eres < 0.) eres = 0.;
1268  return eres;
1269 }
1270 
1271 
1272 
1295 
1296 Double_t KVDetector::GetIncidentEnergy(Int_t Z, Int_t A, Double_t delta_e, enum SolType type)
1297 {
1298  // Overrides KVMaterial::GetIncidentEnergy
1299  // Returns incident energy corresponding to energy loss delta_e in active layer of detector for a given nucleus.
1300  // If delta_e is not given, the current energy loss in the active layer is used.
1301  //
1302  // By default the solution corresponding to the highest incident energy is returned
1303  // This is the solution found for Einc greater than the maximum of the dE(Einc) curve.
1304  // If you want the low energy solution set SolType = KVIonRangeTable::kEmin.
1305  //
1306  // WARNING: calculating the incident energy of a particle using only the dE in a detector
1307  // is ambiguous, as in general (and especially for very heavy ions) the maximum of the dE
1308  // curve occurs for Einc greater than the punch-through energy, therefore it is not always
1309  // true to assume that if the particle does not stop in the detector the required solution
1310  // is that for type=KVIonRangeTable::kEmax. For a range of energies between punch-through
1311  // and dE_max, the required solution is still that for type=KVIonRangeTable::kEmin.
1312  // If the residual energy of the particle is unknown, there is no way to know which is the
1313  // correct solution.
1314  //
1315  // WARNING 2
1316  // If the given energy loss in the active layer is greater than the maximum theoretical dE
1317  // for given Z & A, (dE > GetMaxDeltaE(Z,A)) then we return a NEGATIVE incident energy
1318  // corresponding to the maximum, GetEIncOfMaxDeltaE(Z,A)
1319 
1320  if (Z < 1) return 0.;
1321 
1322  Double_t DE = (delta_e > 0 ? delta_e : GetEnergyLoss());
1323 
1324  // If the given energy loss in the active layer is greater than the maximum theoretical dE
1325  // for given Z & A, (dE > GetMaxDeltaE(Z,A)) then we return a NEGATIVE incident energy
1326  // corresponding to the maximum, GetEIncOfMaxDeltaE(Z,A)
1327  if (DE > GetMaxDeltaE(Z, A)) return -GetEIncOfMaxDeltaE(Z, A);
1328 
1329  TF1* dE = GetELossFunction(Z, A);
1330  Double_t e1, e2;
1331  dE->GetRange(e1, e2);
1332  switch (type) {
1333  case kEmin:
1334  e2 = GetEIncOfMaxDeltaE(Z, A);
1335  break;
1336  case kEmax:
1337  e1 = GetEIncOfMaxDeltaE(Z, A);
1338  break;
1339  }
1340  int status;
1341  Double_t EINC = ProtectedGetX(dE, DE, status, e1, e2);
1342  if (status != 0) {
1343 // Warning("GetIncidentEnergy",
1344 // "In %s : Called for Z=%d A=%d with dE=%f sol_type %d",
1345 // GetName(), Z, A, DE, type);
1346 // Warning("GetIncidentEnergy",
1347 // "Max DeltaE in this case is %f MeV",
1348 // GetMaxDeltaE(Z, A));
1349 // Warning("GetIncidentEnergy",
1350 // "Between Einc limits [%f,%f] no solution found",
1351 // e1,e2);
1352 // Warning("GetIncidentEnergy",
1353 // "Returned value is -Einc for %s dE, %f",
1354 // (status>0 ? "maximum" : "minimum"), EINC);
1355  return -EINC;
1356  }
1357  return EINC;
1358 }
1359 
1360 
1361 
1367 
1368 Double_t KVDetector::GetDeltaEFromERes(Int_t Z, Int_t A, Double_t Eres)
1369 {
1370  // Overrides KVMaterial::GetDeltaEFromERes
1371  //
1372  // Calculate energy loss in active layer of detGetAlignedDetector for nucleus (Z,A)
1373  // having a residual kinetic energy Eres (MeV)
1374 
1375  if (Z < 1 || Eres <= 0.) return 0.;
1376  Double_t Einc = GetIncidentEnergyFromERes(Z, A, Eres);
1377  if (Einc <= 0.) return 0.;
1378  return GetELossFunction(Z, A)->Eval(Einc);
1379 }
1380 
1381 
1382 
1389 
1390 Double_t KVDetector::GetIncidentEnergyFromERes(Int_t Z, Int_t A, Double_t Eres)
1391 {
1392  // Overrides KVMaterial::GetIncidentEnergyFromERes
1393  //
1394  // Calculate incident energy of nucleus from residual energy.
1395  //
1396  // Returns -1 if Eres is out of defined range of values
1397 
1398  if (Z < 1 || Eres <= 0.) return 0.;
1399  //return GetEResFunction(Z, A)->GetX(Eres);
1400  Int_t status;
1401  Double_t einc = KVBase::ProtectedGetX(GetEResFunction(Z, A), Eres, status);
1402  if (status != 0) {
1403  // problem with inversion - value out of defined range of function
1404  return -1;
1405  }
1406  return einc;
1407 }
1408 
1409 
1410 
1414 
1415 Double_t KVDetector::GetSmallestEmaxValid(Int_t Z, Int_t A) const
1416 {
1417  // Returns the smallest maximum energy for which range tables are valid
1418  // for all absorbers in the detector, and given ion (Z,A)
1419 
1420  Double_t maxmin = -1.;
1421  TIter next(&fAbsorbers);
1422  KVMaterial* abs;
1423  while ((abs = (KVMaterial*)next())) {
1424  if (maxmin < 0.) maxmin = abs->GetEmaxValid(Z, A);
1425  else {
1426  if (abs->GetEmaxValid(Z, A) < maxmin) maxmin = abs->GetEmaxValid(Z, A);
1427  }
1428  }
1429  return maxmin;
1430 }
1431 
1432 
1433 
1451 
1452 void KVDetector::ReadDefinitionFromFile(const Char_t* envrc)
1453 {
1454  // Create detector from text file in 'TEnv' format.
1455  //
1456  // Example:
1457  // ========
1458  //
1459  // Layer: Gold
1460  // Gold.Material: Au
1461  // Gold.AreaDensity: 200.*KVUnits::ug
1462  // +Layer: Gas1
1463  // Gas1.Material: C3F8
1464  // Gas1.Thickness: 5.*KVUnits::cm
1465  // Gas1.Pressure: 50.*KVUnits::mbar
1466  // Gas1.Active: yes
1467  // +Layer: Si1
1468  // Si1.Material: Si
1469  // Si1.Thickness: 300.*KVUnits::um
1470 
1471  TEnv fEnvFile(envrc);
1472 
1473  KVString layers(fEnvFile.GetValue("Layer", ""));
1474  layers.Begin(" ");
1475  while (!layers.End()) {
1476  KVString layer = layers.Next();
1477  KVString mat = fEnvFile.GetValue(Form("%s.Material", layer.Data()), "");
1478  KVString tS = fEnvFile.GetValue(Form("%s.Thickness", layer.Data()), "");
1479  KVString pS = fEnvFile.GetValue(Form("%s.Pressure", layer.Data()), "");
1480  KVString dS = fEnvFile.GetValue(Form("%s.AreaDensity", layer.Data()), "");
1481  Double_t thick, dens, press;
1482  thick = dens = press = 0.;
1483  KVMaterial* M = 0;
1484  if (pS != "" && tS != "") {
1485  press = (Double_t)gROOT->ProcessLineFast(Form("%s*1.e+12", pS.Data()));
1486  press /= 1.e+12;
1487  thick = (Double_t)gROOT->ProcessLineFast(Form("%s*1.e+12", tS.Data()));
1488  thick /= 1.e+12;
1489  M = new KVMaterial(mat.Data(), thick, press);
1490  }
1491  else if (tS != "") {
1492  thick = (Double_t)gROOT->ProcessLineFast(Form("%s*1.e+12", tS.Data()));
1493  thick /= 1.e+12;
1494  M = new KVMaterial(mat.Data(), thick);
1495  }
1496  else if (dS != "") {
1497  dens = (Double_t)gROOT->ProcessLineFast(Form("%s*1.e+12", dS.Data()));
1498  dens /= 1.e+12;
1499  M = new KVMaterial(dens, mat.Data());
1500  }
1501  if (M) {
1502  AddAbsorber(M);
1503  if (fEnvFile.GetValue(Form("%s.Active", layer.Data()), kFALSE)) SetActiveLayer(M);
1504  }
1505  }
1506 }
1507 
1508 
1509 
1513 
1514 Double_t KVDetector::GetRange(Int_t Z, Int_t A, Double_t Einc)
1515 {
1516  // WARNING: SAME AS KVDetector::GetLinearRange
1517  // Only linear range in centimetres is calculated for detectors!
1518  return GetLinearRange(Z, A, Einc);
1519 }
1520 
1521 
1522 
1528 
1529 Double_t KVDetector::GetLinearRange(Int_t Z, Int_t A, Double_t Einc)
1530 {
1531  // Returns range of ion in centimetres in this detector,
1532  // taking into account all layers.
1533  // Note that for Einc > punch through energy, this range is no longer correct
1534  // (but still > total thickness of detector).
1535  return GetRangeFunction(Z, A)->Eval(Einc);
1536 }
1537 
1538 
1539 
1543 
1544 Double_t KVDetector::GetPunchThroughEnergy(Int_t Z, Int_t A)
1545 {
1546  // Returns energy (in MeV) necessary for ion (Z,A) to punch through all
1547  // layers of this detector
1548 
1549  if (fSingleLayer) {
1550  // Optimize calculation time for single-layer detector
1551  return fActiveLayer->GetPunchThroughEnergy(Z, A);
1552  }
1553  return GetRangeFunction(Z, A)->GetX(GetTotalThicknessInCM());
1554 }
1555 
1556 
1557 
1558 
1563 
1564 TGraph* KVDetector::DrawPunchThroughEnergyVsZ(Int_t massform)
1565 {
1566  // Creates and fills a TGraph with the punch through energy in MeV vs. Z for the given detector,
1567  // for Z=1-92. The mass of each nucleus is calculated according to the given mass formula
1568  // (see KVNucleus).
1569 
1570  TGraph* punch = new TGraph(92);
1571  punch->SetName(Form("KVDetpunchthrough_%s_mass%d", GetName(), massform));
1572  punch->SetTitle(Form("Simple Punch-through %s (MeV) (mass formula %d)", GetName(), massform));
1573  KVNucleus nuc;
1574  nuc.SetMassFormula(massform);
1575  for (int Z = 1; Z <= 92; Z++) {
1576  nuc.SetZ(Z);
1577  punch->SetPoint(Z - 1, Z, GetPunchThroughEnergy(nuc.GetZ(), nuc.GetA()));
1578  }
1579  return punch;
1580 }
1581 
1582 
1583 
1588 
1589 TGraph* KVDetector::DrawPunchThroughEsurAVsZ(Int_t massform)
1590 {
1591  // Creates and fills a TGraph with the punch through energy in MeV/nucleon vs. Z for the given detector,
1592  // for Z=1-92. The mass of each nucleus is calculated according to the given mass formula
1593  // (see KVNucleus).
1594 
1595  TGraph* punch = new TGraph(92);
1596  punch->SetName(Form("KVDetpunchthroughEsurA_%s_mass%d", GetName(), massform));
1597  punch->SetTitle(Form("Simple Punch-through %s (AMeV) (mass formula %d)", GetName(), massform));
1598  KVNucleus nuc;
1599  nuc.SetMassFormula(massform);
1600  for (int Z = 1; Z <= 92; Z++) {
1601  nuc.SetZ(Z);
1602  punch->SetPoint(Z - 1, Z, GetPunchThroughEnergy(nuc.GetZ(), nuc.GetA()) / nuc.GetA());
1603  }
1604  return punch;
1605 }
1606 
1607 
1608 
1609 
1611 
1612 KVGroup* KVDetector::GetGroup() const
1613 {
1614  return (KVGroup*)GetParentStructure("GROUP");
1615 }
1616 
1617 
1618 
1620 
1621 UInt_t KVDetector::GetGroupNumber()
1622 {
1623  return (GetGroup() ? GetGroup()->GetNumber() : 0);
1624 }
1625 
1626 
1627 
1629 
1630 void KVDetector::AddParentStructure(KVGeoStrucElement* elem)
1631 {
1632  fParentStrucList.Add(elem);
1633 }
1634 
1635 
1636 
1638 
1639 void KVDetector::RemoveParentStructure(KVGeoStrucElement* elem)
1640 {
1641  fParentStrucList.Remove(elem);
1642 }
1643 
1644 
1645 
1649 
1650 KVGeoStrucElement* KVDetector::GetParentStructure(const Char_t* type, const Char_t* name) const
1651 {
1652  // Get parent geometry structure element of given type.
1653  // Give unique name of structure if more than one element of same type is possible.
1654  KVGeoStrucElement* el = 0;
1655  if (strcmp(name, "")) {
1656  KVSeqCollection* strucs = fParentStrucList.GetSubListWithType(type);
1657  el = (KVGeoStrucElement*)strucs->FindObject(name);
1658  delete strucs;
1659  }
1660  else
1661  el = (KVGeoStrucElement*)fParentStrucList.FindObjectByType(type);
1662  return el;
1663 }
1664 
1665 
1666 
1669 
1670 void KVDetector::SetActiveLayerMatrix(const TGeoHMatrix* m)
1671 {
1672  // Set ROOT geometry global matrix transformation to coordinate frame of active layer volume
1673  SetMatrix(m);
1674 }
1675 
1676 
1677 
1680 
1681 void KVDetector::SetActiveLayerShape(TGeoBBox* s)
1682 {
1683  // Set ROOT geometry shape of active layer volume
1684  SetShape(s);
1685 }
1686 
1687 
1688 
1691 
1692 void KVDetector::SetEntranceWindowMatrix(const TGeoHMatrix* m)
1693 {
1694  // Set ROOT geometry global matrix transformation to coordinate frame of entrance window
1695  fEWPosition.SetMatrix(m);
1696 }
1697 
1698 
1699 
1702 
1703 void KVDetector::SetEntranceWindowShape(TGeoBBox* s)
1704 {
1705  // Set ROOT geometry shape of entrance window
1706  fEWPosition.SetShape(s);
1707 }
1708 
1709 
1710 
1721 
1722 void KVDetector::SetThickness(Double_t thick)
1723 {
1724  // Overrides KVMaterial::SetThickness
1725  //
1726  // If ROOT geometry is defined, we modify the DZ thickness of the volume representing
1727  // this detector in accordance with the new thickness.
1728  //
1729  // This is only implemented for single-layer detectors with the following shapes:
1730  // - TGeoBBox (rectangular box)
1731  // - TGeoTube (tube)
1732  // - TGeoArb8 (arbitrary trapezoid with less than 8 vertices standing on two parallel planes perpendicular to Z axis)
1733 
1734  if (ROOTGeo() && fSingleLayer) {
1735  TGeoPhysicalNode* pn = (TGeoPhysicalNode*)gGeoManager->GetListOfPhysicalNodes()->FindObject(GetNode()->GetFullPathToNode());
1736  if (!pn) pn = gGeoManager->MakePhysicalNode(GetNode()->GetFullPathToNode());
1737  TGeoBBox* shape = (TGeoBBox*)pn->GetShape();
1738  TGeoShape* newshape = nullptr;
1739  // bad kludge - is there no better way to clone a shape and change its dZ?
1740  if (shape->IsA() == TGeoBBox::Class()) {
1741  newshape = new TGeoBBox(shape->GetDX(), shape->GetDY(), 0.5 * thick);
1742  }
1743  else if (shape->IsA() == TGeoTube::Class()) {
1744  TGeoTube* oldtube = static_cast<TGeoTube*>(shape);
1745  newshape = new TGeoTube(oldtube->GetRmin(), oldtube->GetRmax(), 0.5 * thick);
1746  }
1747  else if (shape->IsA() == TGeoArb8::Class()) {
1748  TGeoArb8* oldtube = static_cast<TGeoArb8*>(shape);
1749  auto vert = oldtube->GetVertices();
1750  newshape = new TGeoArb8(0.5 * thick, vert);
1751  }
1752  else {
1753  Error("SetThickness", "No implementation for %s (%s)", shape->IsA()->GetName(), GetName());
1754  }
1755  if (newshape) {
1756  pn->Align(nullptr, newshape);
1758  }
1759  }
1760  KVMaterial::SetThickness(thick);
1761 }
1762 
1763 
1764 
1770 
1771 Bool_t KVDetector::HasSameStructureAs(const KVDetector* other) const
1772 {
1773  // Return kTRUE if the two detectors have the same internal structure, i.e.
1774  // - the same number of absorber layers
1775  // - in the same order
1776  // - with the same material & thickness
1777 
1778  int nabs = GetNumberOfAbsorberLayers();
1779  if (other->GetNumberOfAbsorberLayers() != nabs) return kFALSE;
1780  bool same = true;
1781  for (int iabs = 0; iabs < nabs; ++iabs) {
1782  KVMaterial* this_abs = GetAbsorber(iabs);
1783  KVMaterial* other_abs = other->GetAbsorber(iabs);
1784  if (!this_abs->IsType(other_abs->GetType())
1785  || this_abs->GetMass() != other_abs->GetMass()
1786  || this_abs->GetThickness() != other_abs->GetThickness())
1787  same = false;
1788  }
1789  return same;
1790 }
1791 
1792 
1793 
1801 
1802 Bool_t KVDetector::AddDetectorSignalExpression(const KVString& type, const KVString& _expr)
1803 {
1804  // Add a new KVDetectorSignalExpression to this detector
1805  //
1806  // \param[in] type the name/type of the new signal
1807  // \param[in] _expr mathematical expression using any of the known signals of the detector
1808  //
1809  // \note If the expression is not valid, no signal will be created and method returns kFALSE.
1810 
1811  KVDetectorSignalExpression* ds = new KVDetectorSignalExpression(type, _expr, this);
1812  if (!ds->IsValid()) {
1813  delete ds;
1814  ds = nullptr;
1815  }
1816  else
1817  AddDetectorSignal(ds);
1818  return (ds != nullptr);
1819 }
1820 
1821 
1822 
1831 
1832 void KVDetector::SetPressure(Double_t P)
1833 {
1834  // \param[in] P pressure in [Torr]
1835  //
1836  // For a gaseous detector, set/change the pressure of the active gas layer.
1837  //
1838  // For ROOT geometries, we change the medium/material of the corresponding node in the geometry
1839  // in order to reflect the change in pressure (gases are represented by different media/materials
1840  // for each pressure/temperature) so that it will be taken into account for example when filtering simulated data.
1841 
1843  if (ROOTGeo()) {
1844  // find node in geometry
1845  TGeoPhysicalNode* pn = (TGeoPhysicalNode*)gGeoManager->GetListOfPhysicalNodes()->FindObject(GetNode()->GetFullPathToNode());
1846  if (!pn) pn = gGeoManager->MakePhysicalNode(GetNode()->GetFullPathToNode());
1847  // use new medium reflecting change in pressure
1848  pn->GetVolume()->SetMedium(GetActiveLayer()->GetGeoMedium());
1849  }
1850 }
1851 
1852 
1853 
1862 
1863 void KVDetector::SetTemperature(Double_t T)
1864 {
1865  // \param[in] T pressure in [deg. C]
1866  //
1867  // For a gaseous detector, set/change the temperature of the active gas layer.
1868  //
1869  // For ROOT geometries, we change the medium/material of the corresponding node in the geometry
1870  // in order to reflect the change in temperature (gases are represented by different media/materials
1871  // for each pressure/temperature) so that it will be taken into account for example when filtering simulated data.
1872 
1874  if (ROOTGeo()) {
1875  // find node in geometry
1876  TGeoPhysicalNode* pn = (TGeoPhysicalNode*)gGeoManager->GetListOfPhysicalNodes()->FindObject(GetNode()->GetFullPathToNode());
1877  if (!pn) pn = gGeoManager->MakePhysicalNode(GetNode()->GetFullPathToNode());
1878  // use new medium reflecting change in pressure
1879  pn->GetVolume()->SetMedium(GetActiveLayer()->GetGeoMedium());
1880  }
1881 }
1882 
1883 
int Int_t
unsigned int UInt_t
#define SafeDelete(p)
#define e(i)
bool Bool_t
unsigned char UChar_t
char Char_t
float Float_t
constexpr Bool_t kFALSE
double Double_t
constexpr Bool_t kTRUE
const char Option_t
R__EXTERN TEnv * gEnv
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t option
R__EXTERN TGeoManager * gGeoManager
#define gROOT
char * Form(const char *fmt,...)
virtual const Char_t * GetType() const
Definition: KVBase.h:176
virtual Bool_t IsType(const Char_t *typ) const
Definition: KVBase.h:184
Double_t ProtectedGetX(const TF1 *func, Double_t val, int &status, Double_t xmin=0.0, Double_t xmax=0.0) const
Definition: KVBase.cpp:1598
Output signal from detector obtained by calibration.
Bool_t IsAvailableFor(const KVNameValueList &params) const override
Base class for all detector calibrations.
Definition: KVCalibrator.h:99
TString GetInputSignalType() const
Definition: KVCalibrator.h:224
void SetDetector(KVDetector *d)
Definition: KVCalibrator.h:232
TString GetOutputSignalType() const
Definition: KVCalibrator.h:228
Signal output from a mathematical combination of other signals.
Base class for output signal data produced by a detector.
virtual Bool_t IsExpression() const
virtual Bool_t IsRaw() const
virtual Bool_t GetValueNeedsExtraParameters() const
virtual Double_t GetValue(const KVNameValueList &params="") const
virtual void Reset()
Base class describing elements of array geometry.
Group of detectors which can be treated independently of all others in array.
Definition: KVGroup.h:20
Extended TList class which owns its objects by default.
Definition: KVList.h:28
A stack of materials in which successive energy losses of charged particles can be calculated ,...
Description of physical materials used to construct detectors & targets; interface to range tables.
Definition: KVMaterial.h:89
virtual void SetPressure(Double_t)
Definition: KVMaterial.cpp:580
virtual void SetTemperature(Double_t)
Definition: KVMaterial.cpp:646
virtual void SetThickness(Double_t thick)
Definition: KVMaterial.cpp:451
virtual Double_t GetThickness() const
Definition: KVMaterial.cpp:484
void Clear(Option_t *opt="") override
Reset absorber - set stored energy lost by particles in absorber to zero.
virtual Double_t GetERes(Int_t Z, Int_t A, Double_t Einc, Double_t dx=0.)
virtual Double_t GetDeltaE(Int_t Z, Int_t A, Double_t Einc, Double_t dx=0.)
Definition: KVMaterial.cpp:833
virtual Double_t GetLinearRange(Int_t Z, Int_t A, Double_t Einc)
Definition: KVMaterial.cpp:901
Double_t GetMass() const
Definition: KVMaterial.cpp:302
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 * GetStringValue(const Char_t *name) const
Bool_t HasDoubleParameter(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
Int_t GetA() const
Definition: KVNucleus.cpp:792
void SetZ(Int_t z, Char_t mt=-1)
Definition: KVNucleus.cpp:696
void SetMassFormula(UChar_t mt)
Definition: KVNucleus.h:341
Int_t GetZ() const
Return the number of proton / atomic number.
Definition: KVNucleus.cpp:763
TVector3 * GetPInitial() const
Definition: KVParticle.h:729
TVector3 GetMomentum() const
Definition: KVParticle.h:607
KVNameValueList * GetParameters() const
Definition: KVParticle.h:818
Double_t GetEnergy() const
Definition: KVParticle.h:624
void SetE0(TVector3 *e=0)
Definition: KVParticle.h:718
void SetIsDetected()
Definition: KVParticle.h:733
Double_t GetKE() const
Definition: KVParticle.h:617
void SetEnergy(Double_t e)
Definition: KVParticle.h:602
Base class used for handling geometry in a multidetector array.
Definition: KVPosition.h:91
KaliVeda extensions to ROOT collection classes.
TObject * FindObject(const char *name) const override
KVSeqCollection * GetSubListWithType(const Char_t *retvalue) const
Extension of ROOT TString class which allows backwards compatibility with ROOT v3....
Definition: KVString.h:73
KVString Next(Bool_t strip_whitespace=kFALSE) const
Definition: KVString.cpp:695
KVString RNext(Bool_t strip_whitespace=kFALSE) const
Definition: KVString.cpp:841
Handle several calibrations valid for different Z ranges.
virtual const char * GetValue(const char *name, const char *dflt) const
virtual void GetRange(Double_t &xmin, Double_t &xmax) const
Double_t * GetVertices()
static TClass * Class()
virtual Double_t GetDX() const
virtual Double_t GetDY() const
static TClass * Class()
TClass * IsA() const override
void RefreshPhysicalNodes(Bool_t lock=kTRUE)
TGeoPhysicalNode * MakePhysicalNode(const char *path=nullptr)
TObjArray * GetListOfPhysicalNodes()
Bool_t Align(TGeoMatrix *newmat=nullptr, TGeoShape *newshape=nullptr, Bool_t check=kFALSE, Double_t ovlp=0.001)
TGeoVolume * GetVolume(Int_t level=-1) const
TGeoShape * GetShape(Int_t level=-1) const
virtual Double_t GetRmin() const
static TClass * Class()
virtual Double_t GetRmax() const
virtual void SetMedium(TGeoMedium *medium)
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
void SetName(const char *name="") override
void SetTitle(const char *title="") override
TObject * Clone(const char *newname="") const override
const char * GetName() const override
virtual void SetName(const char *name)
TObject * FindObject(const char *name) const override
Longptr_t ExecPlugin(int nargs)
const char * Data() const
RooCmdArg ClassName(const char *name)
RVec< PromoteType< T > > abs(const RVec< T > &v)
Double_t x[n]
gr SetName("gr")
double T(double x)
void Error(const char *location, const char *fmt,...)
void Warning(const char *location, const char *fmt,...)
void init()
constexpr Double_t K()
TArc a
ClassImp(TPyArg)