1 #include "KVDetector.h"
4 #include "KVCalibrator.h"
14 #include "KVCalibratedSignal.h"
15 #include "KVDetectorSignalExpression.h"
16 #include "KVZDependentCalibratedSignal.h"
17 #include "KVMaterialStack.h"
25 Int_t KVDetector::fDetCounter = 0;
31 void KVDetector::init()
34 fCalibrators =
nullptr;
36 fActiveLayer =
nullptr;
37 fIdentP = fUnidentP = 0;
38 fELossF = fEResF = fRangeF =
nullptr;
43 fParentStrucList.SetCleanup();
45 fNode.SetDetector(
this);
47 fDetSignals.SetOwner();
50 fDetSignals.ReplaceObjects();
58 KVDetector::KVDetector()
71 KVDetector::KVDetector(
const Char_t* type,
100 AddAbsorber(
new KVMaterial(gas, thick, pressure, temperature));
114 #if ROOT_VERSION_CODE >= ROOT_VERSION(3,4,0)
117 ((KVDetector&) obj).Copy(*
this);
122 #if ROOT_VERSION_CODE >= ROOT_VERSION(3,4,0)
129 void KVDetector::Copy(
TObject& obj)
const
131 void KVDetector::Copy(
TObject& obj)
138 TIter next(&fAbsorbers);
144 Int_t in_actif = fAbsorbers.IndexOf(fActiveLayer);
145 ((KVDetector&) obj).SetActiveLayer(((KVDetector&)obj).GetAbsorber(in_actif));
153 KVDetector::~KVDetector()
169 void KVDetector::SetMaterial(
const Char_t* type)
175 if (!GetActiveLayer())
178 GetActiveLayer()->SetMaterial(type);
208 if (kvp->
GetKE() <= 0.)
219 std::vector<Double_t> thickness;
224 thickness.reserve(fAbsorbers.GetEntries());
227 TIter next(&fAbsorbers);
229 thickness[i++] =
abs->GetThickness();
231 abs->SetThickness(T);
240 TIter next(&fAbsorbers);
242 abs->SetThickness(thickness[i++]);
251 Double_t eloss_old = GetEnergyLoss();
252 SetEnergyLoss(eloss_old + dE);
279 std::vector<Double_t> thickness;
284 thickness.reserve(fAbsorbers.GetEntries());
287 TIter next(&fAbsorbers);
289 thickness[i++] =
abs->GetThickness();
291 abs->SetThickness(T);
299 TIter next(&fAbsorbers);
301 abs->SetThickness(thickness[i++]);
337 TIter next(&fAbsorbers, kIterBackward);
341 Einc =
abs->GetParticleEIncFromERes(&clone_part, norm);
344 clone_part.SetKE(Einc);
356 void KVDetector::Print(
Option_t* opt)
const
361 if (!strcmp(opt,
"data")) {
362 cout << GetName() <<
" -- ";
365 TIter it(&GetListOfDetectorSignals());
379 if (BelongsToUnidentifiedParticle())
380 cout <<
"(Belongs to an unidentified particle)";
383 else if (!strcmp(opt,
"all")) {
391 TIter next(&fAbsorbers);
393 if (GetActiveLayer() == abs)
394 cout <<
" ### ACTIVE LAYER ### " << endl;
397 if (GetActiveLayer() == abs)
398 cout <<
" #################### " << endl;
401 cout <<
option <<
" --- Detected particles:" << endl;
407 cout <<
ClassName() <<
": " << ((KVDetector*)
this)->
470 Error(
"AddCalibrator",
"Detector:%s has no signal %s, required as input by calibrator %s to provide output %s",
475 fCalibrators =
new KVList();
478 fCalibrators->Add(cal);
482 Warning(
"AddCalibrator",
"%s : output signal not defined for calibrator %s. No output signal created.",
499 if (new_cal_sig) fDetSignals.Add(new_cal_sig);
531 fCalibrators->Remove(old_cal);
532 remove_signal_for_calibrator(old_cal);
535 return AddCalibrator(cal, opts);
565 void KVDetector::Clear(
Option_t* opt)
573 bool option_N = (strncmp(opt,
"N", 1) == 0);
574 bool sim_mode_option_N = option_N && IsSimMode();
578 fIdentP = fUnidentP = 0;
579 ResetBit(kIdentifiedParticle);
580 ResetBit(kUnidentifiedParticle);
582 TIter it(&GetListOfDetectorSignals());
589 if (!sim_mode_option_N) {
592 TIter next(&fAbsorbers);
613 if (!TestBit(kActiveSet))
615 if (fAbsorbers.GetSize() > 1) fSingleLayer =
kFALSE;
627 if (!fAbsorbers.GetEntries()) {
628 Error(
"GetAbsorber",
"No absorbers defined for detector");
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());
646 void KVDetector::RemoveAllAbsorbers()
660 fCalibrators =
nullptr;
661 fParticles =
nullptr;
662 fActiveLayer =
nullptr;
663 ResetBit(kActiveSet);
664 fIdentP = fUnidentP = 0;
665 fELossF = fEResF = fRangeF =
nullptr;
670 fSingleLayer =
kTRUE;
679 void KVDetector::remove_signal_for_calibrator(
KVCalibrator* K)
684 if (
K->GetOutputSignalType() !=
"") {
687 fDetSignals.Remove(ds);
701 void KVDetector::RemoveCalibrators()
710 TIter it(fCalibrators);
712 remove_signal_for_calibrator(K);
714 fCalibrators->Delete();
785 if (e < 0.)
e = GetEnergy();
787 SetEResAfterDetector(-1.);
797 Double_t maxDE = stack.GetMaxDeltaE(z, a);
799 auto inc_angle = stack.GetMinimumIncidentAngleForDEMax(z, a, e);
800 stack.SetIncidenceAngle(inc_angle);
809 return get_corrected_energy(&stack, nuc, e, transmission);
812 Double_t maxDE = GetMaxDeltaE(z, a);
815 auto inc_angle = stack.GetMinimumIncidentAngleForDEMax(z, a, e);
816 stack.SetIncidenceAngle(inc_angle);
824 return get_corrected_energy(&stack, nuc, e, transmission);
826 return get_corrected_energy(
this, nuc, e, transmission);
869 ELOSS = (ELOSS < 0. ? GetEnergy() : ELOSS);
870 if (ELOSS <= 0)
return 0;
877 UInt_t last_positive_difference_z = 1;
879 if (mass_formula > -1)
884 while (zmax > zmin + 1) {
888 difference = GetMaxDeltaE(z, particle.
GetA()) - ELOSS;
890 if (difference < 0.0) {
893 z += (
UInt_t)((zmax - z) / 2 + 0.5);
899 last_positive_difference_z = z;
900 z -= (
UInt_t)((z - zmin) / 2 + 0.5);
904 if (difference < 0) {
907 z = last_positive_difference_z;
935 TIter next(&fAbsorbers);
938 while (fActiveLayer != mat) {
940 e = mat->
GetERes(par[0], par[1], e);
946 return mat->
GetDeltaE(par[0], par[1], e);
977 TIter next(&fAbsorbers);
991 Einc = mat->
GetERes(Z, A, Einc);
1030 TIter next(&fAbsorbers);
1059 KVDetector* KVDetector::MakeDetector(
const Char_t* name,
Float_t thickness)
1077 if (!nom.Contains(
".") && !(ph = LoadPlugin(
"KVDetector", name)))
return nullptr;
1078 if (nom.Contains(
".")) {
1083 if (!(ph = LoadPlugin(
"KVDetector", name))) {
1084 if (!(ph = LoadPlugin(
"KVDetector", det_type))) {
1091 return (KVDetector*) ph->
ExecPlugin(1, thickness);
1099 Double_t KVDetector::GetEntranceWindowSurfaceArea()
1103 if (!GetEntranceWindow().GetShape()->InheritsFrom(
"TGeoArb8")) {
1105 return GetEntranceWindow().GetShape()->GetFacetArea(1);
1108 return GetEntranceWindow().GetSurfaceArea();
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));
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));
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));
1155 fRangeF->SetRange(0., GetSmallestEmaxValid(Z, A));
1156 fRangeF->SetTitle(
Form(
"Range [cm] in detector %s for Z=%d A=%d", GetName(), Z, A));
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));
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));
1198 return GetELossFunction(Z, A)->GetMaximumX();
1214 return GetELossFunction(Z, A)->GetMaximum();
1231 return fActiveLayer->GetDeltaE(Z, A, Einc);
1233 return GetELossFunction(Z, A)->Eval(Einc);
1247 return Einc - GetERes(Z, A, Einc);
1263 if (Einc <= 0.)
return 0.;
1264 Double_t eres = GetEResFunction(Z, A)->Eval(Einc);
1267 if (eres < 0.) eres = 0.;
1320 if (Z < 1)
return 0.;
1322 Double_t DE = (delta_e > 0 ? delta_e : GetEnergyLoss());
1327 if (DE > GetMaxDeltaE(Z, A))
return -GetEIncOfMaxDeltaE(Z, A);
1329 TF1* dE = GetELossFunction(Z, A);
1334 e2 = GetEIncOfMaxDeltaE(Z, A);
1337 e1 = GetEIncOfMaxDeltaE(Z, A);
1341 Double_t EINC = ProtectedGetX(dE, DE, status, e1, e2);
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);
1398 if (Z < 1 || Eres <= 0.)
return 0.;
1421 TIter next(&fAbsorbers);
1424 if (maxmin < 0.) maxmin =
abs->GetEmaxValid(Z, A);
1426 if (
abs->GetEmaxValid(Z, A) < maxmin) maxmin =
abs->GetEmaxValid(Z, A);
1452 void KVDetector::ReadDefinitionFromFile(
const Char_t* envrc)
1471 TEnv fEnvFile(envrc);
1473 KVString layers(fEnvFile.GetValue(
"Layer",
""));
1475 while (!layers.End()) {
1482 thick = dens = press = 0.;
1484 if (pS !=
"" && tS !=
"") {
1491 else if (tS !=
"") {
1496 else if (dS !=
"") {
1503 if (fEnvFile.GetValue(
Form(
"%s.Active", layer.
Data()), kFALSE)) SetActiveLayer(M);
1518 return GetLinearRange(Z, A, Einc);
1535 return GetRangeFunction(Z, A)->Eval(Einc);
1551 return fActiveLayer->GetPunchThroughEnergy(Z, A);
1553 return GetRangeFunction(Z, A)->GetX(GetTotalThicknessInCM());
1564 TGraph* KVDetector::DrawPunchThroughEnergyVsZ(
Int_t massform)
1571 punch->
SetName(
Form(
"KVDetpunchthrough_%s_mass%d", GetName(), massform));
1572 punch->
SetTitle(
Form(
"Simple Punch-through %s (MeV) (mass formula %d)", GetName(), massform));
1575 for (
int Z = 1; Z <= 92; Z++) {
1589 TGraph* KVDetector::DrawPunchThroughEsurAVsZ(
Int_t massform)
1596 punch->
SetName(
Form(
"KVDetpunchthroughEsurA_%s_mass%d", GetName(), massform));
1597 punch->
SetTitle(
Form(
"Simple Punch-through %s (AMeV) (mass formula %d)", GetName(), massform));
1600 for (
int Z = 1; Z <= 92; Z++) {
1612 KVGroup* KVDetector::GetGroup()
const
1614 return (
KVGroup*)GetParentStructure(
"GROUP");
1621 UInt_t KVDetector::GetGroupNumber()
1623 return (GetGroup() ? GetGroup()->GetNumber() : 0);
1632 fParentStrucList.Add(elem);
1641 fParentStrucList.Remove(elem);
1655 if (strcmp(name,
"")) {
1670 void KVDetector::SetActiveLayerMatrix(
const TGeoHMatrix* m)
1681 void KVDetector::SetActiveLayerShape(
TGeoBBox* s)
1692 void KVDetector::SetEntranceWindowMatrix(
const TGeoHMatrix* m)
1695 fEWPosition.SetMatrix(m);
1703 void KVDetector::SetEntranceWindowShape(
TGeoBBox* s)
1706 fEWPosition.SetShape(s);
1722 void KVDetector::SetThickness(
Double_t thick)
1734 if (ROOTGeo() && fSingleLayer) {
1750 newshape =
new TGeoArb8(0.5 * thick, vert);
1753 Error(
"SetThickness",
"No implementation for %s (%s)", shape->
IsA()->
GetName(), GetName());
1756 pn->
Align(
nullptr, newshape);
1771 Bool_t KVDetector::HasSameStructureAs(
const KVDetector* other)
const
1778 int nabs = GetNumberOfAbsorberLayers();
1779 if (other->GetNumberOfAbsorberLayers() != nabs)
return kFALSE;
1781 for (
int iabs = 0; iabs < nabs; ++iabs) {
1783 KVMaterial* other_abs = other->GetAbsorber(iabs);
1817 AddDetectorSignal(ds);
1818 return (ds !=
nullptr);
1832 void KVDetector::SetPressure(
Double_t P)
1863 void KVDetector::SetTemperature(
Double_t T)
winID h TVirtualViewer3D TVirtualGLPainter p
R__EXTERN TGeoManager * gGeoManager
char * Form(const char *fmt,...)
virtual const Char_t * GetType() const
virtual Bool_t IsType(const Char_t *typ) const
Double_t ProtectedGetX(const TF1 *func, Double_t val, int &status, Double_t xmin=0.0, Double_t xmax=0.0) const
Output signal from detector obtained by calibration.
Bool_t IsAvailableFor(const KVNameValueList ¶ms) const override
Base class for all detector calibrations.
TString GetInputSignalType() const
void SetDetector(KVDetector *d)
TString GetOutputSignalType() const
Signal output from a mathematical combination of other signals.
Bool_t IsValid() const override
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 ¶ms="") const
Base class describing elements of array geometry.
Group of detectors which can be treated independently of all others in array.
Extended TList class which owns its objects by default.
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.
virtual void SetPressure(Double_t)
virtual void SetTemperature(Double_t)
virtual void SetThickness(Double_t thick)
virtual Double_t GetThickness() const
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.)
virtual Double_t GetLinearRange(Int_t Z, Int_t A, Double_t Einc)
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.
void SetZ(Int_t z, Char_t mt=-1)
void SetMassFormula(UChar_t mt)
Int_t GetZ() const
Return the number of proton / atomic number.
TVector3 * GetPInitial() const
TVector3 GetMomentum() const
KVNameValueList * GetParameters() const
Double_t GetEnergy() const
void SetE0(TVector3 *e=0)
void SetEnergy(Double_t e)
Base class used for handling geometry in a multidetector array.
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....
KVString Next(Bool_t strip_whitespace=kFALSE) const
KVString RNext(Bool_t strip_whitespace=kFALSE) const
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
virtual Double_t GetDX() const
virtual Double_t GetDY() const
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
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)