kaliveda-dev team mailing list archive
-
kaliveda-dev team
-
Mailing list archive
-
Message #00127
[Merge] lp:~kaliveda-dev/kaliveda/kvedaloss into lp:kaliveda
John Frankland has proposed merging lp:~kaliveda-dev/kaliveda/kvedaloss into lp:kaliveda.
Requested reviews:
Eric Bonnet (ebonnet)
For more details, see:
https://code.launchpad.net/~kaliveda-dev/kaliveda/kvedaloss/+merge/53265
new implementation of VEDALOSS routines (range & energy loss)
separates VEDALOSS from KVMaterial/KVDetector
improved methods for correcting dE for pulse-height defect (silicon detectors) or dead-zones (ChIo etc.)
--
https://code.launchpad.net/~kaliveda-dev/kaliveda/kvedaloss/+merge/53265
Your team KaliVeda Development Team is subscribed to branch lp:kaliveda.
=== modified file 'KVIndra/detectors/KVBIC.cpp'
--- KVIndra/detectors/KVBIC.cpp 2009-11-10 04:41:34 +0000
+++ KVIndra/detectors/KVBIC.cpp 2011-03-14 15:58:00 +0000
@@ -80,24 +80,18 @@
fBomb = bomb;
//build detector
- AddAbsorber(new KVMaterial("Myl", 1.5)); //mylar entry window
+ AddAbsorber(new KVMaterial("Myl", 1.5*KVUnits::um)); //mylar entry window
Float_t e = GetEffectiveEntryThickness();
- KVMaterial *gas = new KVMaterial("CF4", e);
- gas->SetPressure(pressure);
- gas->SetUnits(KVMaterial::kTORR); //set units for pressure to Torr
+ KVMaterial *gas = new KVMaterial("CF4", e*KVUnits::mm, pressure);
AddAbsorber(gas); //gas between two windows
- AddAbsorber(new KVMaterial("Myl", 1.0)); //interior window
- gas = new KVMaterial("CF4", 60);
- gas->SetPressure(pressure);
- gas->SetUnits(KVMaterial::kTORR); //set units for pressure to Torr
+ AddAbsorber(new KVMaterial("Myl", 1.0*KVUnits::um)); //interior window
+ gas = new KVMaterial("CF4", 60*KVUnits::mm,pressure);
AddAbsorber(gas); //main body of gas
SetActiveLayer(gas); //active layer - energy loss is measured
- AddAbsorber(new KVMaterial("Myl", 1.0)); //2nd interor window
- gas = new KVMaterial("CF4", e);
- gas->SetPressure(pressure);
- gas->SetUnits(KVMaterial::kTORR); //set units for pressure to Torr
+ AddAbsorber(new KVMaterial("Myl", 1.0*KVUnits::um)); //2nd interor window
+ gas = new KVMaterial("CF4", e*KVUnits::mm, pressure);
AddAbsorber(gas); //gas between two exit windows
- AddAbsorber(new KVMaterial("Myl", 1.5)); //exit window
+ AddAbsorber(new KVMaterial("Myl", 1.5*KVUnits::um)); //exit window
}
@@ -194,8 +188,8 @@
Bool_t maxDE = kFALSE;
//egas > max possible ?
- if (egas > GetBraggDE(z, a)) {
- egas = GetBraggDE(z, a);
+ if (egas > GetMaxDeltaE(z, a)) {
+ egas = GetMaxDeltaE(z, a);
maxDE = kTRUE;
}
//calculate incident energy
=== modified file 'KVIndra/detectors/KVBIC.h'
--- KVIndra/detectors/KVBIC.h 2009-11-10 04:41:34 +0000
+++ KVIndra/detectors/KVBIC.h 2011-03-14 15:58:00 +0000
@@ -55,6 +55,16 @@
{
return (fLinCal && fLinCal->GetStatus());
};
+ virtual void SetPressure(Double_t P /* torr*/)
+ {
+ // Set pressure of gas in torr
+ GetActiveLayer()->SetPressure(P);
+ };
+ virtual Double_t GetPressure() const /* torr */
+ {
+ // Give pressure of gas in torr
+ return GetActiveLayer()->GetPressure();
+ };
ClassDef(KVBIC, 2) //Blocking ChIo for E416 experiment
};
=== modified file 'KVIndra/detectors/KVChIo.cpp'
--- KVIndra/detectors/KVChIo.cpp 2010-12-03 09:16:52 +0000
+++ KVIndra/detectors/KVChIo.cpp 2011-03-14 15:58:00 +0000
@@ -64,20 +64,19 @@
}
//______________________________________________________________________________
-KVChIo::KVChIo(Float_t pressure, Float_t thick):KVINDRADetector("Myl", 2.5)
+KVChIo::KVChIo(Float_t pressure, Float_t thick):KVINDRADetector("Myl", 2.5*KVUnits::um)
{
- //Make an INDRA ChIo: 2.5micron mylar windows enclosing 'thick' mm of C3F8.
- //By default 'thick'=50mm
- //By default 'pressure'=0mbar
- //The type of these detectors is "CI"
+ // Make an INDRA ChIo: 2.5micron mylar windows enclosing 'thick' cm of C3F8,
+ // give gas pressure in mbar
+ // By default 'thick'=5cm
+ // The type of these detectors is "CI"
//gas layer
- KVMaterial *mat = new KVMaterial("C3F8", thick);
- mat->SetPressure(pressure);
+ KVMaterial *mat = new KVMaterial("C3F8", thick, pressure*KVUnits::mbar);
AddAbsorber(mat);
SetActiveLayer(mat); //gas is the active layer
// mylar for second window
- mat = new KVMaterial("Myl", 2.5);
+ mat = new KVMaterial("Myl", 2.5*KVUnits::um);
AddAbsorber(mat);
SetType("CI");
@@ -181,76 +180,21 @@
//if stopped=kTRUE, we give the correction for a particle which stops in the detector
//(by default we assume the particle continues after the detector)
//
- //If the dE energy loss in the gas is > maximum theoretical dE (GetBraggDE)
- //this calculation is not valid. The mylar energy loss is calculated for a dE
- //corresponding to dE = GetBraggDE(z,a), and then we scale up according
- //to: dE_Mylar = dE_Mylar_Bragg * (dE_measured / dE_Bragg).
- //if dE_measured - dE_Bragg > 2 MeV, a warning is printed.
+ // WARNING: if stopped=kFALSEE, and if the residual energy after the detector
+ // is known (i.e. measured in a detector placed after this one), you should
+ // first call
+ // SetEResAfterDetector(Eres);
+ // before calling this method. Otherwise, especially for heavy ions, the
+ // correction may be false for particles which are just above the punch-through energy.
egas = ((egas < 0.) ? GetEnergy() : egas);
if (egas <= 0.)
return 0.0; //check measured (calibrated) energy in gas is reasonable (>0)
- Bool_t maxDE = kFALSE;
-
- //egas > max possible ?
- Double_t de_measured = 0.;
- if (egas > GetBraggDE(z, a)) {
- de_measured = egas;
- egas = GetBraggDE(z, a);
- maxDE = kTRUE;
-// if(de_measured-egas>2.0)
-// Warning("GetELossMylar","%s: Measured de_ChIo (%f) is greater than maximum for Z=%d (%f)",
-// GetName(), de_measured, (int)z, egas);
- }
- enum KVMaterial::SolType solution = KVMaterial::kEmax;
- if(stopped) solution = KVMaterial::kEmin;
- //calculate incident energy
- Double_t e_inc = GetIncidentEnergy(z, a, egas, solution);
-
- //calculate residual energy
- Double_t e_res = GetERes(z, a, e_inc);
-
- //calculate mylar energy
- Double_t emylar = e_inc - e_res - egas;
-
- emylar = ((emylar < 0.) ? 0. : emylar);
-
- if (maxDE){
- //If the dE energy loss in the gas is > maximum theoretical dE (GetBraggDE)
- //this calculation is not valid. The mylar energy loss is calculated for a dE
- //corresponding to dE = GetBraggDE(z,a), and then we scale up according
- //to: dE_Mylar = dE_Mylar_Bragg * (dE_measured / dE_Bragg)
- emylar *= (de_measured/egas);
- }
-
+ Double_t emylar = GetCorrectedEnergy(z,a,egas,!stopped) - egas;
return emylar;
}
-//__________________________________________________________________________________________________________________________
-
-Double_t KVChIo::GetCorrectedEnergy(UInt_t z, UInt_t a, Double_t egas, Bool_t transmission)
-{
- //Redefinition of KVDetector method.
- //Based on energy loss in gas, calculates correction for mylar windows
- //from energy loss tables. Returns total energy loss in mylar+gas+mylar
- //If argument 'egas' not given, KVChIo::GetEnergy() is used
- //If transmission=kFALSE we give the correction for a particle stopping in the
- //detector (i.e. we calculate the incident energy for a particle with dE<dE_Bragg)
- //By default transmission=kTRUE & we assume the particle continues after the
- //detector.
- //
- //If the dE energy loss in the gas is > maximum theoretical dE (GetBraggDE)
- //this calculation is not valid. The mylar energy loss is calculated for a dE
- //corresponding to dE = GetBraggDE(z,a), and then we scale up according
- //to: dE_Mylar = dE_Mylar_Bragg * (dE_measured / dE_Bragg)
-
- egas = ((egas < 0.) ? GetEnergy() : egas);
- if( egas <=0 ) return 0;
- Double_t emyl = GetELossMylar(z, a, egas, !transmission);
- return (egas + emyl);
-}
-
//_______________________________________________________________________________
Double_t KVChIo::GetCalibratedEnergy()
=== modified file 'KVIndra/detectors/KVChIo.h'
--- KVIndra/detectors/KVChIo.h 2010-12-03 09:16:52 +0000
+++ KVIndra/detectors/KVChIo.h 2011-03-14 15:58:00 +0000
@@ -41,8 +41,7 @@
public:
KVChIo();
- KVChIo(Float_t pressure, Float_t thick = 50.0);
-// KVChIo(Float_t thick = 50.0);
+ KVChIo(Float_t pressure, Float_t thick = 5.0*KVUnits::cm);
virtual ~ KVChIo();
Double_t GetVoltsFromCanalPG(Double_t chan = 0.0);
@@ -77,11 +76,21 @@
};
Double_t GetELossMylar(UInt_t z, UInt_t a, Double_t egas = -1.0, Bool_t stopped=kFALSE);
- Double_t GetCorrectedEnergy(UInt_t z, UInt_t a, Double_t egas = -1.0, Bool_t transmission=kTRUE);
inline Bool_t IsCalibrated() const;
virtual Short_t GetCalcACQParam(KVACQParam*,Double_t) const;
+ virtual void SetPressure(Double_t P /* mbar */)
+ {
+ // Set pressure of gas in mbar
+ GetActiveLayer()->SetPressure(P*KVUnits::mbar);
+ };
+ virtual Double_t GetPressure() const /* mbar */
+ {
+ // Give pressure of gas in mbar
+ return GetActiveLayer()->GetPressure()/KVUnits::mbar;
+ };
+
ClassDef(KVChIo, 4) //The ionisation chamber detectors (ChIo) of the INDRA array
};
=== modified file 'KVIndra/detectors/KVChIo_e475s.h'
--- KVIndra/detectors/KVChIo_e475s.h 2010-12-03 09:16:52 +0000
+++ KVIndra/detectors/KVChIo_e475s.h 2011-03-14 15:58:00 +0000
@@ -25,7 +25,7 @@
public:
KVChIo_e475s();
- KVChIo_e475s(Float_t pressure, Float_t thick=50.0);
+ KVChIo_e475s(Float_t pressure, Float_t thick=50.0*KVUnits::mm);
virtual ~KVChIo_e475s(){};
virtual void SetCalibrators(){};
=== modified file 'KVIndra/detectors/KVPhoswich.cpp'
--- KVIndra/detectors/KVPhoswich.cpp 2009-11-02 08:25:48 +0000
+++ KVIndra/detectors/KVPhoswich.cpp 2011-03-14 15:58:00 +0000
@@ -60,10 +60,6 @@
KVPhoswich::~KVPhoswich()
{
- //If gIndra exists, remove from list of phoswich
-// if (gIndra) {
-// gIndra->GetListOfPhoswich()->Remove(this);
-// }
}
=== modified file 'KVIndra/detectors/KVSiB.cpp'
--- KVIndra/detectors/KVSiB.cpp 2009-11-10 04:41:34 +0000
+++ KVIndra/detectors/KVSiB.cpp 2011-03-14 15:58:00 +0000
@@ -32,7 +32,7 @@
KVSiB::KVSiB(Float_t thickness):KVSilicon(thickness),fLinCal(0)
{
- //Make a 'thickness' um thick silicon detector with type="SIB"
+ // Make a 'thickness' um thick silicon detector with type="SIB"
SetType("SIB");
}
=== modified file 'KVIndra/detectors/KVSilicon.cpp'
--- KVIndra/detectors/KVSilicon.cpp 2010-12-03 09:16:52 +0000
+++ KVIndra/detectors/KVSilicon.cpp 2011-03-14 15:58:00 +0000
@@ -39,7 +39,7 @@
//
//Used to describe Silicon detectors of the INDRA array.
//In order to create a detector, use the KVSilicon::KVSilicon(Float_t thick)
-//constructor with "thick" the thickness in microns.
+//constructor with "thick" the thickness in centimetres
//
//Type of detector : "SI"
@@ -70,10 +70,10 @@
}
//______________________________________________________________________________
-KVSilicon::KVSilicon(Float_t thick):KVINDRADetector("Si", thick)
+KVSilicon::KVSilicon(Float_t thick):KVINDRADetector("Si", thick*KVUnits::um)
{
- // constructor for silicon detector
- //Type of detector: "SI"
+ // constructor for silicon detector, thickness in microns
+ // Type of detector: "SI"
SetType("SI");
init();
}
@@ -235,30 +235,6 @@
return fPHD->Compute(Einc);
}
-//____________________________________________________________________________________________
-
-Double_t KVSilicon::GetCorrectedEnergy(UInt_t z, UInt_t a, Double_t e, Bool_t trn)
-{
- //Returns total energy lost by particle in silicon detector corrected for PHD.
- //
- //If "e" (measured/apparent energy loss in detector) not given, current value
- //measured in detector is used. If PHD for detector has not been set, no correction
- //is performed
- //See GetPHD() and class KVPulseHeightDefect.
- //
- //transmission=kTRUE (default): particle does not stop in the detector
- //transmission=kFALSE: particle stops in the detector
-
- if (e < 0.) e = GetEnergy();
- if( e <= 0 ) return 0;
- // calculate incident energy from measured energy loss in detector
- enum KVMaterial::SolType solution = KVMaterial::kEmax;
- if(!trn) solution = KVMaterial::kEmin;
- Double_t EINC = GetIncidentEnergy(z, a, e, solution);
- //incident energy - residual energy = total real energy loss
- return (EINC - GetERes(z, a, EINC));
-}
-
//__________________________________________________________________________________________
Double_t KVSilicon::GetCalibratedEnergy()
{
@@ -414,63 +390,6 @@
//______________________________________________________________________________
-void KVSilicon::SetELossParams(Int_t Z, Int_t A)
-{
- //Initialise energy loss coefficients for this detector and a given incident nucleus (Z,A)
- //We redefine the KVDetector::SetELossParams method in order to include the
- //pulse height defect (if defined) for the silicon detector in the calculation of
- //the energy loss in the active (silicon) layer.
-
- //PHD not defined ? ignore
- if(!fPHD || !fPHD->GetStatus()){
- KVDetector::SetELossParams(Z,A);
- return;
- }
-
- //do we need to set up the ELoss function ?
- //only if it has never been done for PHD before
- Int_t npar_siphd = 19+6; //number of params for eloss function = 19 params for energy loss + 6 for PHD (Z,a_1,a_2,b_1,b_2,Zmin)
-
- if( npar_loss != npar_siphd ){
- npar_loss = npar_siphd;
- if( par_loss ) delete [] par_loss; //delete previous parameter array
- par_loss = new Double_t[npar_loss];
- //find/create function
- //search in list of functions for one corresponding to this detector
- //the name of the required function is ELoss_SiPHD
- ELoss =
- (TF1 *) gROOT->GetListOfFunctions()->FindObject("ELoss_SiPHD");
- if (!ELoss)
- ELoss = new TF1("ELoss_SiPHD", ELossSiPHD, 0.1, 5000., npar_loss);
- }
-
- //fill parameter array
- ((KVMaterial*)fAbsorbers->At(0))->GetELossParams(Z, A, par_loss);
- par_loss[19] = Z;
- for (register int i = 0; i < 5; i++) par_loss[i+20] = fPHD->GetParameter(i);
-
- //set parameters of energy loss function
- ELoss->SetParameters(par_loss);
-}
-
-//_____________________________________________________________________________________//
-
-Double_t ELossSiPHD(Double_t * x, Double_t * par)
-{
- //Calculates measured/apparent energy loss in silicon detector, taking into account PHD
- //Parameters par[0] to par [18] are energy loss parameters for Silicon
- //Parameter par[19] is Z of nucleus
- //Parameters par[20] to par[24] are Moulton PHD parameters (see KVPulseHeightDefect
- //source file for function & parameter definitions).
- //Argument x[0] is incident energy in MeV
-
- //measured/apparent dE = real dE - PHD
- return (ELossSaclay(x, par) - PHDMoulton(x, &par[19]));
-}
-
-
-//______________________________________________________________________________
-
Short_t KVSilicon::GetCalcACQParam(KVACQParam* ACQ,Double_t ECalc) const
{
// Calculates & returns value of given acquisition parameter corresponding to
@@ -484,6 +403,26 @@
return -1;
}
+//______________________________________________________________________________
+
+TF1* KVSilicon::GetELossFunction(Int_t Z, Int_t A)
+{
+ // Overrides KVDetector::GetELossFunction
+ // If the pulse height deficit (PHD) has been set for this detector,
+ // we return an energy loss function which takes into account the PHD,
+ // i.e. for an incident energy E we calculate dEphd(E,Z,A) = dE(E,Z,A) - PHD(E',Z)
+ // (where E' is the energy just before the active layer, in case there are
+ // dead zones before it)
+ // If no PHD is set, we return the usual KVDetector::GetELossFunction
+ // which calculates dE(E,Z,A)
+
+ if(fPHD && fPHD->GetStatus()) {
+ fELossF = fPHD->GetELossFunction(Z,A);
+ fELossF->SetRange(0., GetSmallestEmaxValid(Z,A));
+ }
+ return KVDetector::GetELossFunction(Z,A);
+}
+
//__________________________________________________________________________________________
ClassImp(KVSi75)
@@ -541,7 +480,7 @@
//Default ctor
// first layer (active) : 2mm silicon (nominal)
// second layer (dead) : 40um silicon (nominal)
- AddAbsorber( new KVMaterial("Si", 40.0) );
+ AddAbsorber( new KVMaterial("Si", 40.0*KVUnits::um) );
SetType("SILI");
SetLabel("SILI");
}
=== modified file 'KVIndra/detectors/KVSilicon.h'
--- KVIndra/detectors/KVSilicon.h 2010-12-03 09:16:52 +0000
+++ KVIndra/detectors/KVSilicon.h 2011-03-14 15:58:00 +0000
@@ -29,8 +29,6 @@
class KVChIo;
class KVDBParameterSet;
-Double_t ELossSiPHD(Double_t * x, Double_t * par);
-
class KVSilicon:public KVINDRADetector {
protected:
@@ -52,7 +50,7 @@
public:
KVSilicon();
- KVSilicon(Float_t thick);
+ KVSilicon(Float_t thick /* um */);
virtual ~ KVSilicon();
Double_t GetVoltsFromCanalPG(Double_t chan = 0.0);
@@ -90,8 +88,6 @@
};
Double_t GetPHD(Double_t Einc, UInt_t Z);
- Double_t GetCorrectedEnergy(UInt_t z, UInt_t a, Double_t e = -1., Bool_t transmission=kTRUE);
- virtual void SetELossParams(Int_t Z, Int_t A);
inline Bool_t IsCalibrated() const;
@@ -100,7 +96,19 @@
void SetZminPHD(Int_t zmin) { fZminPHD = zmin; };
Int_t GetZminPHD() { return fZminPHD; };
virtual Short_t GetCalcACQParam(KVACQParam*,Double_t) const;
+ virtual TF1* GetELossFunction(Int_t Z, Int_t A);
+ virtual void SetThickness(Double_t thick /* um */)
+ {
+ // Sets thickness of active layer in microns
+ GetActiveLayer()->SetThickness(thick*KVUnits::um);
+ };
+ virtual Double_t GetThickness() const /* um */
+ {
+ // Returns thickness of active layer in microns
+ return GetActiveLayer()->GetThickness()/KVUnits::um;
+ };
+
ClassDef(KVSilicon, 7) //INDRA forward-rings silicon detector
};
=== modified file 'KVIndra/particles/KVINDRAReconNuc.cpp'
--- KVIndra/particles/KVINDRAReconNuc.cpp 2011-01-21 16:40:26 +0000
+++ KVIndra/particles/KVINDRAReconNuc.cpp 2011-03-14 15:58:00 +0000
@@ -846,6 +846,7 @@
// therefore we have to estimate the silicon energy for this particle using the CsI energy
if(!fPileup && fCoherent && GetSi()->IsCalibrated()){
// all is apparently well
+ GetSi()->SetEResAfterDetector(fECsI);
fESi = GetSi()->GetCorrectedEnergy(GetZ(),GetA());
if(fESi<=0) {
SetECode(kECode15);// bad - no Si energy
@@ -856,6 +857,7 @@
{
Double_t e0 = GetSi()->GetDeltaEFromERes(GetZ(),GetA(),fECsI);
// calculated energy: negative
+ GetSi()->SetEResAfterDetector(fECsI);
fESi = -GetSi()->GetCorrectedEnergy(GetZ(),GetA(),e0);
SetECode(kECode2);
}
@@ -865,12 +867,14 @@
// if fUseFullChIoEnergyForCalib = kFALSE, we have to estimate the ChIo energy for this particle
if(fUseFullChIoEnergyForCalib && GetChIo()->IsCalibrated()){
// all is apparently well
+ GetChIo()->SetEResAfterDetector(TMath::Abs(fESi)+fECsI);
fEChIo = GetChIo()->GetCorrectedEnergy(GetZ(),GetA());
}
else
{
- Double_t e0 = GetChIo()->GetDeltaEFromERes(GetZ(),GetA(),fECsI+fESi);
+ Double_t e0 = GetChIo()->GetDeltaEFromERes(GetZ(),GetA(),TMath::Abs(fESi)+fECsI);
// calculated energy: negative
+ GetChIo()->SetEResAfterDetector(TMath::Abs(fESi)+fECsI);
fEChIo = -GetChIo()->GetCorrectedEnergy(GetZ(),GetA(),e0);
SetECode(kECode2);
}
=== modified file 'KVMultiDet/KVMultiDetLinkDef.h'
--- KVMultiDet/KVMultiDetLinkDef.h 2011-03-14 08:44:04 +0000
+++ KVMultiDet/KVMultiDetLinkDef.h 2011-03-14 15:58:00 +0000
@@ -17,6 +17,7 @@
#pragma link C++ enum KVBase::EKaliVedaBits;
#pragma link C++ function SearchFile(const Char_t*, TString&, int);
#pragma link C++ namespace KVTGIDFunctions;
+#pragma link C++ namespace KVUnits;
#pragma link C++ nestedclass;
#pragma link C++ nestedtypedef;
#pragma extra_include "Rtypes.h";
@@ -126,7 +127,7 @@
#pragma link C++ class KVLVContainer;
#pragma link C++ class KVLVEntry;
#pragma link C++ class KVLVColumnData;
-#pragma link C++ class KVMaterial-;
+#pragma link C++ class KVMaterial+;
#pragma link C++ class KVMemoryChunk+;
#pragma link C++ class KVMemoryPool+;
#pragma link C++ class KVModule+;
@@ -253,6 +254,7 @@
#pragma link C++ class SRBFile_t+;
#pragma link C++ class SRBDataRepository+;
#pragma link C++ class SRBAvailableRunsFile+;
+<<<<<<< TREE
#pragma link C++ class KVNuclData+;
#pragma link C++ class KVLifeTime+;
#pragma link C++ class KVAbundance+;
@@ -264,4 +266,9 @@
#pragma link C++ class KVAbundanceTable+;
#pragma link C++ class KVMassExcessTable+;
#pragma link C++ class KVNDTManager+;
+=======
+#pragma link C++ class KVIonRangeTable+;
+#pragma link C++ class KVedaLossMaterial+;
+#pragma link C++ class KVedaLoss+;
+>>>>>>> MERGE-SOURCE
#endif
=== modified file 'KVMultiDet/Makefile'
--- KVMultiDet/Makefile 2009-12-08 08:32:41 +0000
+++ KVMultiDet/Makefile 2011-03-14 15:58:00 +0000
@@ -3,7 +3,7 @@
PROJ_NAME = KVMultiDet
-MODULES = analysis base calibration daq_cec data_management db detectors events geometry gui identification idmaps idtelescopes particles \
+MODULES = analysis base stopping calibration daq_cec data_management db detectors events geometry gui identification idmaps idtelescopes particles \
vg_base vg_charge vg_energy vg_multiplicity vg_shape vg_momentum trieur
EXTRA_HEADERS = base/KVError.h base/KVParameter.h
=== modified file 'KVMultiDet/calibration/KVPulseHeightDefect.cpp'
--- KVMultiDet/calibration/KVPulseHeightDefect.cpp 2009-01-22 15:03:32 +0000
+++ KVMultiDet/calibration/KVPulseHeightDefect.cpp 2011-03-14 15:58:00 +0000
@@ -37,7 +37,7 @@
// --> END_HTML
////////////////////////////////////////////////////////////////////////////////
-Double_t PHDMoulton(Double_t * x, Double_t * par)
+Double_t KVPulseHeightDefect::PHDMoulton(Double_t * x, Double_t * par)
{
//Returns Moulton PHD for given E and Z:
//
@@ -51,21 +51,14 @@
//
// x[0] = E (MeV)
// par[0] = Z
- // par[1] = a_1
- // par[2] = a_2
- // par[3] = b_1
- // par[4] = b_2
- // par[5] = Zmin
- if( par[0] <= par[5] ) return 0;
- Double_t Z = par[0];
- Double_t a = par[1]*Z*Z/1000. + par[2];
- Double_t b = par[3]*100./Z + par[4];
+ Int_t Z = par[0];
+ if( Z <= fZmin ) return 0;
+ Double_t a = a_1*Z*Z/1000. + a_2;
+ Double_t b = b_1*100./Z + b_2;
return (TMath::Power(10,b)*TMath::Power(x[0],a));
}
-TF1 KVPulseHeightDefect::fMoulton("fMoulton", PHDMoulton, 0., 5000., 6);
-
//___________________________________________________________________//
void KVPulseHeightDefect::init()
@@ -73,6 +66,9 @@
//default initialisations
SetType("Pulse Height Defect");
SetZ(1);
+ fMoulton = fDeltaEphd = 0;
+ a_1=a_2=b_1=b_2=0.;
+ fZmin=0;
}
KVPulseHeightDefect::KVPulseHeightDefect():KVCalibrator(5)
@@ -109,19 +105,7 @@
// b(Z) = b_1*(100/Z) + b_2
// E = incident energy of particle
- const_cast<KVPulseHeightDefect*>(this)->SetPars();
- return fMoulton.Eval(energy);
-}
-
-//___________________________________________________________________________
-
-void KVPulseHeightDefect::SetPars()
-{
- // Set parameters of TF1 prior to calculation (PRIVATE)
- Double_t par[6];
- for (register int i = 0; i < 5; i++) par[i+1] = GetParameter(i);
- par[0] = (Double_t) GetZ();
- fMoulton.SetParameters(par);
+ return const_cast<KVPulseHeightDefect*>(this)->GetMoultonPHDFunction(GetZ())->Eval(energy);
}
//___________________________________________________________________________
@@ -133,6 +117,62 @@
return Compute(energy);
}
+TF1* KVPulseHeightDefect::GetELossFunction(Int_t Z, Int_t A)
+{
+ // Return pointer to TF1 giving energy loss in active layer of detector minus
+ // the pulse height defect for a given nucleus (Z,A).
+
+ if(!fDeltaEphd){
+ fDeltaEphd = new TF1(Form("KVPulseHeightDefect:%s:ELossActive", GetDetector()->GetName()),
+ this, &KVPulseHeightDefect::ELossActive, 0., 1.e+04, 2, "KVPulseHeightDefect", "ELossActive");
+ fDeltaEphd->SetNpx( gEnv->GetValue("KVPulseHeightDefect.EnergyLoss.Npx", 20) );
+ }
+ fDeltaEphd->SetParameters((Double_t)Z, (Double_t)A);
+ GetMoultonPHDFunction(Z);
+ return fDeltaEphd;
+}
+
+Double_t KVPulseHeightDefect::ELossActive(Double_t *x, Double_t *par)
+{
+ // Calculate energy lost in active layer of detector minus the Moulton PHD
+ // x[0] = incident energy
+ // par[0] = Z
+ // par[1] = A
+
+ Double_t e = x[0];
+ TIter next(GetDetector()->GetListOfAbsorbers()); KVMaterial* mat;
+ while( (mat = (KVMaterial*)next()) != GetDetector()->GetActiveLayer() ){
+ // calculate energy losses in absorbers before active layer
+ e = mat->GetERes(par[0], par[1], e); //residual energy after layer
+ if (e < KVDETECTOR_MINIMUM_E)
+ return 0.; // return 0 if particle stops in layers before active layer
+ }
+ // calculate energy loss in active layer
+ Double_t dE = mat->GetDeltaE(par[0], par[1], e);
+ // calculate Moulton PHD corresponding to incident energy in active layer
+ Double_t phd = PHDMoulton(&e, par);
+
+ return dE - phd;
+}
+
+ TF1* KVPulseHeightDefect::GetMoultonPHDFunction(Int_t Z)
+ {
+ // Create TF1* fMoulton if not already done
+
+ if(!fMoulton) {
+ fMoulton = new TF1(Form("MoultonPHD:%s", GetDetector()->GetName()),
+ this, &KVPulseHeightDefect::PHDMoulton, 0., 1.e+04, 1, "KVPulseHeightDefect", "PHDMoulton");
+ a_1 = GetParameter(0);
+ a_2 = GetParameter(1);
+ b_1 = GetParameter(2);
+ b_2 = GetParameter(3);
+ fZmin = (Int_t)GetParameter(4);
+ fMoulton->SetNpx(500);
+ }
+ fMoulton->SetParameter(0,Z);
+ return fMoulton;
+ }
+
//___________________________________________________________________________
Double_t KVPulseHeightDefect::Invert(Double_t energy)
@@ -141,7 +181,7 @@
//(set using SetZ() method), this method inverts the Moulton formula
//in order to find the incident energy of the particle.
- const_cast<KVPulseHeightDefect*>(this)->SetPars();
- Double_t xmin, xmax; fMoulton.GetRange(xmin,xmax);
- return fMoulton.GetX(energy,xmin,xmax);
+ const_cast<KVPulseHeightDefect*>(this)->GetMoultonPHDFunction(GetZ());
+ Double_t xmin, xmax; fMoulton->GetRange(xmin,xmax);
+ return fMoulton->GetX(energy,xmin,xmax);
}
=== modified file 'KVMultiDet/calibration/KVPulseHeightDefect.h'
--- KVMultiDet/calibration/KVPulseHeightDefect.h 2009-01-22 15:03:32 +0000
+++ KVMultiDet/calibration/KVPulseHeightDefect.h 2011-03-14 15:58:00 +0000
@@ -12,14 +12,14 @@
#include <KVCalibrator.h>
-Double_t PHDMoulton(Double_t * x, Double_t * par);
class KVPulseHeightDefect : public KVCalibrator
{
- static TF1 fMoulton; //Moulton formula for PHD = f(E,Z)
- UInt_t fZ; //!Z of nucleus to be calibrated
-
- void SetPars();
+ TF1* fMoulton; //!Moulton formula for PHD = f(E,Z)
+ TF1* fDeltaEphd; //!deltaE calculated including PHD
+ Int_t fZ; //!Z of nucleus to be calibrated
+ Int_t fZmin; //!minimum Z for which PHD is considered
+ Double_t a_1, a_2, b_1, b_2; //! parameters of Moulton formula
public:
@@ -33,12 +33,16 @@
virtual Double_t operator() (Double_t );
virtual Double_t Invert(Double_t);
- void SetZ(UInt_t z) {
+ void SetZ(Int_t z) {
fZ = z;
};
- UInt_t GetZ() const {
+ Int_t GetZ() const {
return fZ;
};
+ Double_t PHDMoulton(Double_t * x, Double_t * par);
+ TF1* GetMoultonPHDFunction(Int_t Z);
+ Double_t ELossActive(Double_t *x, Double_t *par);
+ TF1* GetELossFunction(Int_t Z, Int_t A);
ClassDef(KVPulseHeightDefect,1)//Silicon PHD described by Moulton formula
};
=== modified file 'KVMultiDet/data/kvloss.data'
--- KVMultiDet/data/kvloss.data 2010-05-18 08:51:06 +0000
+++ KVMultiDet/data/kvloss.data 2011-03-14 15:58:00 +0000
@@ -1,14 +1,14 @@
//$Id: kvloss.data,v 1.4 2007/06/14 12:29:37 franklan Exp $
//Modified version of loss1.vedasac for KVMaterial classes
//First line gives
-// type-name, name, units flag, density, Z, atomic mass, thickness of gas (mm),
-// molecular weight, temperature [last 3 only used for gases].
-// units flag is for internal purposes and determines how "thickness" is calculated.
+// type-name, name, state, density (g/cm**3), Z, atomic mass,
+// molecular weight, temperature [last 2 only used for gases].
+// 'state' can be 'solid' or 'gas'
//Lines beginning "Z" give limits of validity for energy loss calculation
//as a function of Z and E/A.
-+ Si Silicon 1 2.33 14.0 28.0855 0. 0. 0.
++ Si Silicon solid 2.33 14.0 28.0855 0. 0.
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1391E+1 0.1432E+1 0.1042E+0 -.1066E-1 0.5099E-4 0.4765E-5
-.1073E+1 0.8721E+0 -.8156E-1 0.1146E-1 -.8426E-3 0.2541E-4
2. 4. 0.1500E+1 0.1309E+1 0.1525E+0 -.1717E-1 -.1733E-3 0.7908E-4
@@ -210,9 +210,9 @@
100. 252. 0.1634E+1 0.5772E+0 -.1199E-1 0.1884E-1 0.2472E-2 -.5296E-3
-.2370E+1 0.7501E+0 0.8252E+0 -.2933E+0 0.3903E-1 -.1851E-2
-+ Myl Mylar 1 1.395 4.545 8.735 0. 0. 0.
++ Myl Mylar solid 1.395 4.545 8.735 0. 0.
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.9435E+0 0.1577E+1 0.8671E-1 -.1634E-1 0.1885E-2 -.1336E-3
-.6187E+0 0.7095E+0 -.4930E-1 0.8294E-2 -.7041E-3 0.2359E-4
2. 4. 0.1044E+1 0.1454E+1 0.1395E+0 -.2414E-1 0.1769E-2 -.6188E-4
@@ -414,9 +414,9 @@
100. 252. 0.1182E+1 0.5533E+0 0.1721E-1 0.1818E-1 0.1377E-2 -.3958E-3
-.2066E+1 0.1569E+1 0.3647E+0 -.2140E+0 0.3561E-1 -.2003E-2
-+ NE102 Plastic 3 1.032 3.5 6.509 0. 0. 0.
++ NE102 Plastic solid 1.032 3.5 6.509 0. 0.
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.8330E+0 0.1605E+1 0.8125E-1 -.1629E-1 0.2028E-2 -.1501E-3
-.5293E+0 0.6794E+0 -.4110E-1 0.7086E-2 -.6185E-3 0.2139E-4
2. 4. 0.9288E+0 0.1486E+1 0.1340E+0 -.2498E-1 0.2146E-2 -.9344E-4
@@ -618,9 +618,9 @@
100. 252. 0.1102E+1 0.5582E+0 0.1074E-1 0.1767E-1 0.2284E-2 -.5073E-3
-.1922E+1 0.1616E+1 0.3091E+0 -.1967E+0 0.3324E-1 -.1881E-2
-+ Ni Nickel 2 8.902 28.0 58.6934 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ Ni Nickel solid 8.902 28.0 58.6934 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1794E+1 0.1345E+1 0.9702E-1 -.4056E-2 -.8964E-3 0.4412E-4
-.1508E+1 0.9607E+0 -.8070E-1 0.9007E-2 -.5604E-3 0.1557E-4
2. 4. 0.1892E+1 0.1234E+1 0.1429E+0 -.1142E-1 -.8038E-3 0.9291E-4
@@ -822,9 +822,9 @@
100. 252. 0.2010E+1 0.5539E+0 -.9370E-2 0.2132E-1 0.4241E-3 -.2647E-3
-.2367E+1 -.3687E+0 0.1447E+1 -.4287E+0 0.5236E-1 -.2345E-2
-+ C3F8 C3F8 0 0 8.182 17.093 50.0 188.0202 19.0
++ C3F8 Octofluoropropane gas 0 8.182 17.093 188.0202 19.0
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1187E+1 0.1512E+1 0.9917E-1 -.1521E-1 0.1205E-2 -.7625E-4
-.8428E+0 0.7909E+0 -.6996E-1 0.1100E-1 -.8766E-3 0.2780E-4
2. 4. 0.1297E+1 0.1385E+1 0.1502E+0 -.2191E-1 0.8933E-3 0.9256E-5
@@ -1026,9 +1026,9 @@
100. 252. 0.1387E+1 0.5537E+0 0.6666E-1 0.2016E-1 -.4188E-2 0.2614E-3
-.2498E+1 0.1574E+1 0.3905E+0 -.2295E+0 0.3905E-1 -.2237E-2
-+ C Carbon 2 1.90 6.0 12.011 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ C Carbon solid 1.90 6.0 12.011 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1010E+1 0.1559E+1 0.9390E-1 -.1724E-1 0.1840E-2 -.1249E-3
-.6781E+0 0.7369E+0 -.5892E-1 0.9897E-2 -.8293E-3 0.2731E-4
2. 4. 0.1117E+1 0.1431E+1 0.1464E+0 -.2438E-1 0.1560E-2 -.3813E-4
@@ -1230,9 +1230,9 @@
100. 252. 0.1263E+1 0.5824E+0 0.7019E-2 0.1550E-1 0.2413E-2 -.4678E-3
-.2102E+1 0.1482E+1 0.3437E+0 -.1808E+0 0.2776E-1 -.1451E-2
-+ Ag Silver 2 10.5 47.0 107.8682 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ Ag Silver solid 10.5 47.0 107.8682 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.2170E+1 0.1307E+1 0.8207E-1 -.6085E-3 -.7715E-3 -.1175E-5
-.1877E+1 0.9851E+0 -.6665E-1 0.5782E-2 -.3041E-3 0.8360E-5
2. 4. 0.2253E+1 0.1206E+1 0.1275E+0 -.8997E-2 -.4412E-3 0.3295E-4
@@ -1434,9 +1434,9 @@
100. 252. 0.2256E+1 0.5653E+0 0.1612E-2 0.1933E-1 -.8577E-3 -.5125E-4
-.2326E+1 -.8741E+0 0.1598E+1 -.4352E+0 0.5037E-1 -.2167E-2
-+ Sn Tin 2 5.75 50.0 118.710 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ Sn Tin solid 5.75 50.0 118.710 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.2238E+1 0.1304E+1 0.7985E-1 -.3463E-3 -.7042E-3 -.1061E-4
-.1941E+1 0.9876E+0 -.6487E-1 0.5446E-2 -.2817E-3 0.7822E-5
2. 4. 0.2319E+1 0.1205E+1 0.1252E+0 -.8871E-2 -.3341E-3 0.1999E-4
@@ -1639,9 +1639,9 @@
-.2595E+1 -.5102E+0 0.1384E+1 -.3793E+0 0.4376E-1 -.1875E-2
-+ CsI Cesium_Iodide 3 4.51 54. 129.905 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ CsI CesiumIodide solid 4.51 54. 129.905 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.2294E+1 0.1302E+1 0.7701E-1 -.8533E-4 -.6047E-3 -.2335E-4
-.1990E+1 0.9862E+0 -.6225E-1 0.5043E-2 -.2581E-3 0.7331E-5
2. 4. 0.2372E+1 0.1204E+1 0.1223E+0 -.8750E-2 -.2054E-3 0.5590E-5
@@ -1844,9 +1844,9 @@
-.2908E+1 -.2595E-1 0.1068E+1 -.2878E+0 0.3175E-1 -.1291E-2
-+ Au Gold 2 19.3 79.0 196.96655 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ Au Gold solid 19.3 79.0 196.96655 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.2536E+1 0.1300E+1 0.6188E-1 0.3559E-3 0.8626E-4 -.9808E-4
-.2178E+1 0.9665E+0 -.5087E-1 0.3869E-2 -.2256E-3 0.7655E-5
2. 4. 0.2601E+1 0.1208E+1 0.1075E+0 -.8978E-2 0.6010E-3 -.7457E-4
@@ -2049,9 +2049,9 @@
-.9877E+0 -.3494E+1 0.2679E+1 -.6177E+0 0.6352E-1 -.2468E-2
-+ U Uranium 2 18.95 92. 238.0289 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ U Uranium solid 18.95 92. 238.0289 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.2652E+1 0.1304E+1 0.5572E-1 0.1675E-3 0.4154E-3 -.1292E-3
-.2263E+1 0.9587E+0 -.4833E-1 0.3900E-2 -.2553E-3 0.9003E-5
2. 4. 0.2711E+1 0.1214E+1 0.1016E+0 -.9389E-2 0.9604E-3 -.1066E-3
@@ -2254,9 +2254,9 @@
-.2286E+0 -.4604E+1 0.3112E+1 -.6884E+0 0.6837E-1 -.2567E-2
-+ Air Air 0 0 7.212 14.428 50.0 28.848 19.0
++ Air Air gas 0 7.212 14.428 28.848 19.0
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1086E+1 0.1532E+1 0.9708E-1 -.1632E-1 0.1539E-2 -.1001E-3
-.7506E+0 0.7647E+0 -.6505E-1 0.1061E-1 -.8693E-3 0.2808E-4
2. 4. 0.1194E+1 0.1403E+1 0.1490E+0 -.2301E-1 0.1179E-2 -.1159E-4
@@ -2459,9 +2459,9 @@
-.2494E+1 0.1599E+1 0.4181E+0 -.2492E+0 0.4320E-1 -.2521E-2
-+ Nb Nobium 2 8.57 41.0 92.9064 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ Nb Nobium solid 8.57 41.0 92.9064 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.2082E+1 0.1314E+1 0.8674E-1 -.1315E-2 -.8874E-3 0.1774E-4
-.1795E+1 0.9853E+0 -.7145E-1 0.6646E-2 -.3636E-3 0.9829E-5
2. 4. 0.2169E+1 0.1211E+1 0.1321E+0 -.9419E-2 -.5992E-3 0.5259E-4
@@ -2664,9 +2664,9 @@
-.2307E+1 -.8474E+0 0.1626E+1 -.4498E+0 0.5264E-1 -.2283E-2
-+ Ta Tantalum 2 16.654 73. 180.9479 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ Ta Tantalum solid 16.654 73. 180.9479 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.2488E+1 0.1299E+1 0.6507E-1 0.3753E-3 -.7426E-4 -.8208E-4
-.2144E+1 0.9715E+0 -.5281E-1 0.3983E-2 -.2205E-3 0.7233E-5
2. 4. 0.2556E+1 0.1206E+1 0.1106E+0 -.8829E-2 0.4173E-3 -.5739E-4
@@ -2869,9 +2869,9 @@
-.1138E+1 -.3268E+1 0.2621E+1 -.6159E+0 0.6429E-1 -.2531E-2
-+ Al Aluminium 2 2.6989 13.0 26.9815 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ Al Aluminium solid 2.6989 13.0 26.9815 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1389E+1 0.1443E+1 0.1040E+0 -.1138E-1 0.2061E-3 -.5113E-5
-.1063E+1 0.8660E+0 -.8158E-1 0.1164E-1 -.8629E-3 0.2610E-4
2. 4. 0.1499E+1 0.1319E+1 0.1527E+0 -.1785E-1 -.4496E-4 0.7172E-4
@@ -3074,9 +3074,9 @@
-.2350E+1 0.7300E+0 0.8520E+0 -.3065E+0 0.4146E-1 -.2001E-2
-+ KCl KCl 2 1.987 18.0 37.27565 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ KCl KCl solid 1.987 18.0 37.27565 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1549E+1 0.1397E+1 0.1033E+0 -.8198E-2 -.4010E-3 0.3012E-4
-.1242E+1 0.9134E+0 -.8462E-1 0.1105E-1 -.7693E-3 0.2246E-4
2. 4. 0.1656E+1 0.1278E+1 0.1506E+0 -.1491E-1 -.5298E-3 0.9645E-4
@@ -3279,9 +3279,9 @@
-.2593E+1 0.8169E+0 0.7883E+0 -.2800E+0 0.3712E-1 -.1760E-2
-+ CF4 CF4 0 0 8.40 17.601 50.0 88.046 19.0
++ CF4 Tetrafluoromethane gas 0 8.40 17.601 88.046 19.0
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1201E+1 0.1508E+1 0.9964E-1 -.1504E-1 0.1150E-2 -.7201E-4
-.8563E+0 0.7956E+0 -.7093E-1 0.1110E-1 -.8806E-3 0.2784E-4
2. 4. 0.1310E+1 0.1381E+1 0.1505E+0 -.2171E-1 0.8374E-3 0.1327E-4
@@ -3484,9 +3484,9 @@
-.2530E+1 0.1568E+1 0.4070E+0 -.2358E+0 0.3998E-1 -.2285E-2
-+ Ca Calcium 2 1.550 20.0 40.078 0. 0. 0.
++ Ca Calcium solid 1.550 20.0 40.078 0. 0.
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1570E+1 0.1384E+1 0.1024E+0 -.7159E-2 -.5674E-3 0.3821E-4
-.1271E+1 0.9212E+0 -.8349E-1 0.1058E-1 -.7226E-3 0.2094E-4
2. 4. 0.1675E+1 0.1267E+1 0.1492E+0 -.1404E-1 -.6303E-3 0.9995E-4
@@ -3689,9 +3689,9 @@
-.2509E+1 0.4779E+0 0.9968E+0 -.3303E+0 0.4250E-1 -.1972E-2
-+ Ge Germanium 2 5.323 32.0 72.64 0. 0. 0.
++ Ge Germanium solid 5.323 32.0 72.64 0. 0.
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1946E+1 0.1333E+1 0.9387E-1 -.2968E-2 -.9443E-3 0.3906E-4
-.1661E+1 0.9806E+0 -.7940E-1 0.8308E-2 -.4920E-3 0.1336E-4
2. 4. 0.2040E+1 0.1224E+1 0.1392E+0 -.1052E-1 -.7384E-3 0.7290E-4
@@ -3894,9 +3894,9 @@
-.2308E+1 -.8158E+0 0.1715E+1 -.4976E+0 0.6075E-1 -.2739E-2
-+ Cu Copper 2 8.96 29.0 63.546 0. 0. 0.
++ Cu Copper solid 8.96 29.0 63.546 0. 0.
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1857E+1 0.1342E+1 0.9627E-1 -.3763E-2 -.9160E-3 0.4353E-4
-.1571E+1 0.9701E+0 -.8112E-1 0.8887E-2 -.5444E-3 0.1497E-4
2. 4. 0.1954E+1 0.1231E+1 0.1425E+0 -.1116E-1 -.8653E-3 0.9969E-4
@@ -4099,9 +4099,9 @@
-.2375E+1 -.4707E+0 0.1519E+1 -.4513E+0 0.5569E-1 -.2528E-2
-+ Ti Titanium 2 4.54 22.0 47.867 0. 0. 0.
++ Ti Titanium solid 4.54 22.0 47.867 0. 0.
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1703E+1 0.1372E+1 0.1013E+0 -.6260E-2 -.6836E-3 0.4217E-4
-.1404E+1 0.9461E+0 -.8558E-1 0.1043E-1 -.6891E-3 0.1953E-4
2. 4. 0.1806E+1 0.1256E+1 0.1478E+0 -.1324E-1 -.7105E-3 0.1002E-3
@@ -4304,9 +4304,9 @@
-.2437E+1 -.1083E-1 0.1313E+1 -.4180E+0 0.5378E-1 -.2521E-2
-+ Bi Bismuth 2 9.747 83.0 208.98038 0. 0. 0.
++ Bi Bismuth solid 9.747 83.0 208.98038 0. 0.
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.2572E+1 0.1301E+1 0.5992E-1 0.2980E-3 0.1936E-3 -.1084E-3
-.2204E+1 0.9642E+0 -.5002E-1 0.3867E-2 -.2341E-3 0.8064E-5
2. 4. 0.2666E+1 0.1199E+1 0.1000E+0 -.6415E-2 0.4247E-3 -.8060E-4
@@ -4509,9 +4509,9 @@
-.2321E+1 -.1423E+1 0.1639E+1 -.3809E+0 0.3806E-1 -.1413E-2
-+ V Vanadium 2 6.11 23. 50.9415 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ V Vanadium solid 6.11 23. 50.9415 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1743E+1 0.1367E+1 0.1007E+0 -.5818E-2 -.7426E-3 0.4433E-4
-.1446E+1 0.9535E+0 -.8556E-1 0.1024E-1 -.6671E-3 0.1876E-4
2. 4. 0.1846E+1 0.1252E+1 0.1470E+0 -.1285E-1 -.7460E-3 0.1000E-3
@@ -4714,9 +4714,9 @@
-.2479E+1 0.2175E-1 0.1277E+1 -.4073E+0 0.5263E-1 -.2483E-2
-+ C4H10 Isobutane 0 0 2.429 4.152 50.0 58.0 19.0
++ C4H10 Isobutane gas 0 2.429 4.152 58.12 19.0
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.6534E+0 0.1643E+1 0.7193E-1 -.1573E-1 0.2192E-2 -.1704E-3
-.3984E+0 0.6425E+0 -.3041E-1 0.5452E-2 -.4990E-3 0.1815E-4
2. 4. 0.7417E+0 0.1529E+1 0.1249E+0 -.2546E-1 0.2586E-2 -.1347E-3
@@ -4919,9 +4919,9 @@
-.1660E+1 0.1866E+1 -.2700E-1 -.9656E-1 0.2255E-1 -.1550E-2
-+ Pb Lead 2 11.35 82. 207.2 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ Pb Lead solid 11.35 82. 207.2 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.2569E+1 0.1301E+1 0.6035E-1 0.3209E-3 0.1681E-3 -.1060E-3
-.2203E+1 0.9651E+0 -.5021E-1 0.3860E-2 -.2315E-3 0.7947E-5
2. 4. 0.2669E+1 0.1197E+1 0.9961E-1 -.5993E-2 0.3577E-3 -.7784E-4
@@ -5124,9 +5124,9 @@
-.1936E+1 -.2071E+1 0.1994E+1 -.4702E+0 0.4870E-1 -.1902E-2
-+ PbS Lead_sulphide 2 7.5 49. 120. 0. 0. 0.
-Z = 1,2 0.1 < E/A < 500 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
++ PbS LeadSulphide solid 7.5 49. 120. 0. 0.
+Z = 1,2 0.1 < E/A < 400 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.2332E+1 0.1340E+1 0.6956E-1 -.3771E-2 0.5808E-3 -.1136E-3
-.1977E+1 0.9829E+0 -.7385E-1 0.8338E-2 -.5732E-3 0.1743E-4
2. 4. 0.2427E+1 0.1226E+1 0.1142E+0 -.9551E-2 0.3095E-3 -.3971E-4
@@ -5329,9 +5329,9 @@
-.2333E+1 -.1236E+1 0.1689E+1 -.4269E+0 0.4632E-1 -.1876E-2
-+ Mg Magnesium 2 1.738 12.0 24.3050 0. 0. 0.
++ Mg Magnesium solid 1.738 12.0 24.3050 0. 0.
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1327E+1 0.1454E+1 0.1036E+0 -.1213E-1 0.3751E-3 -.1623E-4
-.1001E+1 0.8495E+0 -.7926E-1 0.1156E-1 -.8722E-3 0.2670E-4
2. 4. 0.1446E+1 0.1320E+1 0.1566E+0 -.1906E-1 0.8074E-4 0.6822E-4
@@ -5534,9 +5534,9 @@
-.2382E+1 0.8949E+0 0.7874E+0 -.3016E+0 0.4242E-1 -.2118E-2
-+Li Lithium 1 0.534 3.0 6.941 0. 0. 0.
++Li Lithium solid 0.534 3.0 6.941 0. 0.
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.9314E+0 0.1634E+1 0.8194E-1 -.1796E-1 0.2213E-2 -.1475E-3
-.5856E+0 0.6761E+0 -.4253E-1 0.7258E-2 -.6097E-3 0.2010E-4
2. 4. 0.1030E+1 0.1514E+1 0.1360E+0 -.2708E-1 0.2363E-2 -.8904E-4
@@ -5738,9 +5738,9 @@
100. 252. 0.1217E+1 0.5142E+0 0.1427E-1 0.2357E-1 0.1572E-2 -.5516E-3
-.2171E+1 0.1448E+1 0.5935E+0 -.3137E+0 0.5230E-1 -.2971E-2
-+Zn Zinc 1 7.133 30.0 65.390 0. 0. 0.
++Zn Zinc solid 7.133 30.0 65.390 0. 0.
Z = 1,2 0.1 < E/A < 400 MeV
-Z = 3,100 0.1 < E/A < 150 MeV
+Z = 3,100 0.1 < E/A < 250 MeV
1. 1. 0.1871E+1 0.1339E+1 0.9546E-1 -.3501E-2 -.9226E-3 0.4174E-4
-.1585E+1 0.9714E+0 -.8023E-1 0.8677E-2 -.5270E-3 0.1447E-4
2. 4. 0.1972E+1 0.1223E+1 0.1439E+0 -.1129E-1 -.8221E-3 0.9275E-4
=== modified file 'KVMultiDet/detectors/KVDetector.cpp'
--- KVMultiDet/detectors/KVDetector.cpp 2011-02-18 14:02:41 +0000
+++ KVMultiDet/detectors/KVDetector.cpp 2011-03-14 15:58:00 +0000
@@ -108,11 +108,11 @@
fIDTelAlign->SetCleanup(kTRUE);
fIDTele4Ident=0;
fIdentP = fUnidentP = 0;
- npar_loss = npar_res = 0;
- par_loss = par_res = 0;
fTotThickness = 0.;
fDepthInTelescope = 0.;
fFiredMask.Set("0");
+ fELossF = fEResF = 0;
+ fEResforEinc = -1.;
}
KVDetector::KVDetector()
@@ -125,7 +125,7 @@
KVDetector::KVDetector(const Char_t * type,
const Float_t thick):KVMaterial()
{
- //Create a new detector of a given material and thickness (default value = 0.0)
+ // Create a new detector of a given material and thickness in centimetres (default value = 0.0)
init();
SetType("DET");
@@ -176,22 +176,22 @@
SafeDelete(fParticles);
delete fAbsorbers;
SafeDelete(fACQParams);
+ SafeDelete(fELossF);
+ SafeDelete(fEResF);
init();
fActiveLayer = -1;
fIDTelAlign->Clear();
SafeDelete(fIDTelAlign);
- if(par_loss) delete [] par_loss;
- if(par_res) delete [] par_res;
SafeDelete(fIDTele4Ident);
-
}
//________________________________________________________________
void KVDetector::SetMaterial(const Char_t * type)
{
- //Set material of active layer.
- //If no absorbers have been added to the detector, create and add
- //one (active layer by default)
+ // Set material of active layer.
+ // If no absorbers have been added to the detector, create and add
+ // one (active layer by default)
+
if (!GetActiveLayer())
AddAbsorber(new KVMaterial(type));
else
@@ -546,6 +546,7 @@
while ((mat = (KVMaterial *) next())) {
mat->Clear();
}
+ fEResforEinc=-1.;
}
//______________________________________________________________________________
@@ -602,7 +603,7 @@
if (fcon == KVD_RECPRC_CNXN)
T->AddDetector(this, KVD_NORECPRC_CNXN);
} else {
- Warning("AddToTelescope", KVDETECTOR_ADD_TO_UNKNOWN_TELESCOPE);
+ Warning("AddToTelescope", "Pointer to telescope is null");
}
}
@@ -808,32 +809,55 @@
Double_t KVDetector::GetCorrectedEnergy(UInt_t z, UInt_t a, Double_t e, Bool_t transmission)
{
- //This function should be redefined in specific detector child classes in order
- //to return the total energy loss in the detector for a given nucleus
- //including inactive absorber layers, and any particle-dependent correction to
- //the 'raw' energy calibration.
+ // Returns the total energy loss in the detector for a given nucleus
+ // including inactive absorber layers.
+ // e = energy loss in active layer (if not given, we use current value)
+ // transmission = kTRUE (default): the particle is assumed to emerge with
+ // a non-zero residual energy Eres after the detector.
+ // = kFALSE: the particle is assumed to stop in the detector.
//
- //If not redefined, this just returns the same as GetEnergy()
+ // WARNING: if transmission=kTRUE, and if the residual energy after the detector
+ // is known (i.e. measured in a detector placed after this one), you should
+ // first call
+ // SetEResAfterDetector(Eres);
+ // before calling this method. Otherwise, especially for heavy ions, the
+ // correction may be false for particles which are just above the punch-through energy.
- if (e > 0.)
- return e;
- return GetEnergy();
+ if (e < 0.) e = GetEnergy();
+ if( e <= 0 ) { SetEResAfterDetector(-1.); return 0; }
+
+ enum KVIonRangeTable::SolType solution = KVIonRangeTable::kEmax;
+ if(!transmission) solution = KVIonRangeTable::kEmin;
+
+ Double_t EINC, ERES = GetEResAfterDetector();
+ if(transmission && ERES>0.){
+ // if residual energy is known we use it to calculate EINC.
+ // if EINC < max of dE curve, we change solution
+ EINC = GetIncidentEnergyFromERes(z, a, ERES);
+ if(EINC < GetEIncOfMaxDeltaE(z,a)) solution = KVIonRangeTable::kEmin;
+ // we could keep the EINC value calculated using ERES, but then
+ // the corrected dE of this detector would not depend on the
+ // measured dE !
+ }
+ EINC = GetIncidentEnergy(z, a, e, solution);
+ ERES = GetERes(z,a,EINC);
+
+ SetEResAfterDetector(-1.);
+ //incident energy - residual energy = total real energy loss
+ return (EINC - ERES);
}
//______________________________________________________________________________//
-UInt_t KVDetector::FindZmin(Double_t ELOSS, Char_t mass_formula)
+Int_t KVDetector::FindZmin(Double_t ELOSS, Char_t mass_formula)
{
//For particles which stop in the first stage of an identification telescope,
//we can at least estimate a minimum Z value based on the energy lost in this
//detector.
//
- //This is based on the KVMaterial::GetBraggDE method, giving the maximum
- //energy loss in the detector for a given nucleus. For e.g. single-layer silicon
- //detectors this energy corresponds to the threshold or punch-through energy
- //(including correction for pulse-height deficit, if known), whereas
- //for e.g. gas detectors with mylar windows this energy is slightly higher than that
- //corresponding to the limit of the particle stopping in the detector.
+ //This is based on the KVMaterial::GetMaxDeltaE method, giving the maximum
+ //energy loss in the active layer of the detector for a given nucleus (A,Z).
+ //
//The "Zmin" is the Z of the particle which gives a maximum energy loss just greater
//than that measured in the detector. Particles with Z<Zmin could not lose as much
//energy and so are excluded.
@@ -864,7 +888,7 @@
particle.SetZ(z);
- difference = GetBraggDE(z,particle.GetA()) - ELOSS;
+ difference = GetMaxDeltaE(z,particle.GetA()) - ELOSS;
//if difference < 0 the z is too small
if (difference < 0.0) {
@@ -890,133 +914,56 @@
//_____________________________________________________________________________________//
-Double_t ELossActive(Double_t * x, Double_t * par)
+Double_t KVDetector::ELossActive(Double_t * x, Double_t * par)
{
- //Calculates energy loss in active layer of detector, taking into account preceding layers
- //First parameter is number of active layer, first layer is numbered 1
- //Following parameters (by sets of 19) are parameters for energy loss in different layers
- //Argument x[0] is incident energy in MeV
-
- Int_t n_active = (Int_t) par[0];
- if (n_active < 2)
- return ELossSaclay(x, &par[1]);
- else {
- Int_t par_ind = 1;
- Double_t e = x[0];
- for (Int_t layer = 1; layer < n_active; layer++) {
-
- e = EResSaclay(&e, &par[par_ind]); //residual energy after layer
+ // Calculates energy loss (in MeV) in active layer of detector, taking into account preceding layers
+ //
+ // Arguments are:
+ // x[0] is incident energy in MeV
+ // Parameters are:
+ // par[0] Z of ion
+ // par[1] A of ion
+
+ Double_t e = x[0];
+ TIter next(fAbsorbers); KVMaterial* mat;
+ if(fActiveLayer>0){
+ // calculate energy losses in absorbers before active layer
+ for (Int_t layer = 0; layer < fActiveLayer; layer++) {
+
+ mat = (KVMaterial*)next();
+ e = mat->GetERes(par[0], par[1], e); //residual energy after layer
if (e < KVDETECTOR_MINIMUM_E)
return 0.; // return 0 if particle stops in layers before active layer
- par_ind += 19;
}
- //calculate energy loss in active layer
- return ELossSaclay(&e, &par[par_ind]);
}
- return 0.;
+ mat = (KVMaterial*)next();
+ //calculate energy loss in active layer
+ return mat->GetDeltaE(par[0], par[1], e);
}
//_____________________________________________________________________________________//
-Double_t EResDet(Double_t * x, Double_t * par)
-{
- //Calculates residual energy of particle after traversing all layers of detector.
- //Returned value is 0 if particle stops in one of the layers of the detector.
- //First parameter par[0]=number of layers of detector
- //Following parameters (by sets of 19) are parameters for energy loss in different layers
- //Argument x[0] is incident energy in MeV
-
- Int_t n_layers = (Int_t) par[0];
- if (n_layers < 2)
- return EResSaclay(x, &par[1]);
- else {
- Int_t par_ind = 1;
- Double_t e = x[0];
- for (Int_t layer = 1; layer <= n_layers; layer++) {
-
- e = EResSaclay(&e, &par[par_ind]); //residual energy after layer
- if (e < KVDETECTOR_MINIMUM_E)
- return 0.; // return 0 if particle stops in one of layers in detector
- par_ind += 19;
-
- }
- return e;
- }
- return 0.;
-}
-
-//______________________________________________________________________________________________//
-
-void KVDetector::SetELossParams(Int_t Z, Int_t A)
-{
- //Initialise energy loss coefficients for this detector and a given incident nucleus (Z,A)
-
- //do we need to set up the ELoss function ?
- //only if it has never been done or if the index of the active layer has changed...
- if( !npar_loss || (npar_loss != 1 + ((Int_t) fActiveLayer + 1) * 19) ){
- //number of params for eloss function = 1 (index of active layer) + 19 params for each layer up to & including active layer
- npar_loss = 1 + ((Int_t) fActiveLayer + 1) * 19;
- if( par_loss ) delete [] par_loss; //delete previous parameter array
- par_loss = new Double_t[npar_loss];
- par_loss[0] = (Int_t) fActiveLayer + 1;
- //find/create function
- TString name_eloss; name_eloss.Form("ELoss_nact%d", ((Int_t) fActiveLayer + 1));
- //search in list of functions for one corresponding to this detector
- //the name of the required function is ELoss_nactx with x = number of active layer (1, 2, etc.)
- ELoss =
- (TF1 *) gROOT->GetListOfFunctions()->FindObject(name_eloss.Data());
- if (!ELoss)
- ELoss = new TF1(name_eloss.Data(), ELossActive, 0.1, 5000., npar_loss);
- }
-
- //loop over layers to fill parameter arrays
- TIter next_abs(fAbsorbers); int i=0; KVMaterial*abs;
- while( (abs = (KVMaterial*)next_abs()) ){
- if (i <= (int) fActiveLayer)
- abs->GetELossParams(Z, A, (par_loss + 1 + 19 * i));
- else break;
- i++;
- }
-
- //set parameters of energy loss function
- ELoss->SetParameters(par_loss);
-}
-
-//______________________________________________________________________________________________//
-
-void KVDetector::SetEResParams(Int_t Z, Int_t A)
-{
- //Initialise coefficients of residual energy function for this detector and a given incident nucleus (Z,A)
-
- //do we need to set up the ERes function ?
- //only if it has never been done or if the number of layers has changed...
- if( !npar_res || (npar_res != 1 + fAbsorbers->GetSize() * 19) ){
- //number of params for eres function = 1 (number of layers) + 19 params for each layer in detector
- npar_res = 1 + fAbsorbers->GetSize() * 19;
- if( par_res ) delete [] par_res; //delete previous parameter array
- par_res = new Double_t[npar_res];
- par_res[0] = fAbsorbers->GetSize();
- //find/create function
- TString name_eres; name_eres.Form("ERes_nlay%d", fAbsorbers->GetSize());
- //search in list of functions for one corresponding to this material
- //the name of the required function is ERes_nlayx with x= number of absorbers in detector
- ERes =
- (TF1 *) gROOT->GetListOfFunctions()->FindObject(name_eres.Data());
- if (!ERes)
- ERes =
- new TF1(name_eres.Data(), EResDet, 0.1, 5000., npar_res);
- }
-
- //loop over layers to fill parameter arrays
- TIter next_abs(fAbsorbers); int i=0; KVMaterial*abs;
- while( (abs = (KVMaterial*)next_abs()) ){
- abs->GetELossParams(Z, A, (par_res + 1 + 19 * i));
- i++;
- }
-
- //set parameters of residual energy function
- ERes->SetParameters(par_res);
+Double_t KVDetector::EResDet(Double_t * x, Double_t * par)
+{
+ // Calculates residual energy (in MeV) of particle after traversing all layers of detector.
+ // Returned value is -1000 if particle stops in one of the layers of the detector.
+ //
+ // Arguments are:
+ // x[0] is incident energy in MeV
+ // Parameters are:
+ // par[0] Z of ion
+ // par[1] A of ion
+
+ Double_t e = x[0];
+ TIter next(fAbsorbers); KVMaterial* mat;
+ while( (mat = (KVMaterial*)next()) ){
+ Double_t eres = mat->GetERes(par[0], par[1], e); //residual energy after layer
+ if (eres <= 0.)
+ return -1000.; // return -1000 if particle stops in layers before active layer
+ e = eres;
+ }
+ return e;
}
//____________________________________________________________________________________
@@ -1146,7 +1093,7 @@
while( (abs = (KVMaterial*)next()) ){
// get medium for absorber
med = abs->GetGeoMedium();
- Double_t thick = abs->GetThicknessInCM();
+ Double_t thick = abs->GetThickness();
GetVerticesInOwnFrame(coords, fDepthInTelescope+depth_in_det, thick);
for(register int i=0;i<8;i++){
vertices[2*i] = coords[i].X();
@@ -1299,4 +1246,210 @@
return surf;
}
+
+TF1* KVDetector::GetEResFunction(Int_t Z, Int_t A)
+{
+ // Return pointer toTF1 giving residual energy after detector as function of incident energy,
+ // for a given nucleus (Z,A).
+ // The TF1::fNpx parameter is taken from environment variable KVDetector.ResidualEnergy.Npx
+
+ if(!fEResF){
+ fEResF = new TF1(Form("KVDetector:%s:ERes", GetName()), this, &KVDetector::EResDet,
+ 0., 1.e+04, 2, "KVDetector", "EResDet");
+ fEResF->SetNpx( gEnv->GetValue("KVDetector.ResidualEnergy.Npx", 20) );
+ }
+ fEResF->SetParameters((Double_t)Z, (Double_t)A);
+ fEResF->SetRange(0., GetSmallestEmaxValid(Z,A));
+
+ return fEResF;
+}
+
+TF1* KVDetector::GetELossFunction(Int_t Z, Int_t A)
+{
+ // Return pointer to TF1 giving energy loss in active layer of detector as function of incident energy,
+ // for a given nucleus (Z,A).
+ // The TF1::fNpx parameter is taken from environment variable KVDetector.EnergyLoss.Npx
+
+ if(!fELossF){
+ fELossF = new TF1(Form("KVDetector:%s:ELossActive", GetName()), this, &KVDetector::ELossActive,
+ 0., 1.e+04, 2, "KVDetector", "ELossActive");
+ fELossF->SetNpx( gEnv->GetValue("KVDetector.EnergyLoss.Npx", 20));
+ }
+ fELossF->SetParameters((Double_t)Z, (Double_t)A);
+ fELossF->SetRange(0., GetSmallestEmaxValid(Z,A));
+ return fELossF;
+}
+
+Double_t KVDetector::GetEIncOfMaxDeltaE(Int_t Z, Int_t A)
+{
+ // Overrides KVMaterial::GetEIncOfMaxDeltaE
+ // Returns incident energy corresponding to maximum energy loss in the
+ // active layer of the detector, for a given nucleus.
+
+ return GetELossFunction(Z,A)->GetMaximumX();
+}
+
+Double_t KVDetector::GetDeltaE(Int_t Z, Int_t A, Double_t Einc)
+{
+ // Overrides KVMaterial::GetDeltaE
+ // Returns energy loss of given nucleus in the active layer of the detector.
+
+ return GetELossFunction(Z,A)->Eval(Einc);
+}
+
+Double_t KVDetector::GetTotalDeltaE(Int_t Z, Int_t A, Double_t Einc)
+{
+ // Returns calculated total energy loss of ion in ALL layers of the detector.
+ // This is just (Einc - GetERes(Z,A,Einc))
+
+ return Einc-GetERes(Z,A,Einc);
+}
+
+Double_t KVDetector::GetERes(Int_t Z, Int_t A, Double_t Einc)
+{
+ // Overrides KVMaterial::GetERes
+ // Returns residual energy of given nucleus after the detector.
+
+ Double_t eres = GetEResFunction(Z,A)->Eval(Einc);
+ // Eres function returns -1000 when particle stops in detector,
+ // in order for function inversion (GetEIncFromEres) to work
+ if(eres<0.) eres = 0.;
+ return eres;
+}
+
+Double_t KVDetector::GetIncidentEnergy(Int_t Z, Int_t A, Double_t delta_e, enum KVIonRangeTable::SolType type)
+{
+ // Overrides KVMaterial::GetIncidentEnergy
+ // Returns incident energy corresponding to energy loss delta_e in active layer of detector for a given nucleus.
+ // If delta_e is not given, the current energy loss in the active layer is used.
+ //
+ // By default the solution corresponding to the highest incident energy is returned
+ // This is the solution found for Einc greater than the maximum of the dE(Einc) curve.
+ // If you want the low energy solution set SolType = KVIonRangeTable::kEmin.
+ //
+ // WARNING: calculating the incident energy of a particle using only the dE in a detector
+ // is ambiguous, as in general (and especially for very heavy ions) the maximum of the dE
+ // curve occurs for Einc greater than the punch-through energy, therefore it is not always
+ // true to assume that if the particle does not stop in the detector the required solution
+ // is that for type=KVIonRangeTable::kEmax. For a range of energies between punch-through
+ // and dE_max, the required solution is still that for type=KVIonRangeTable::kEmin.
+ // If the residual energy of the particle is unknown, there is no way to know which is the
+ // correct solution.
+ //
+ // If the given energy loss in the active layer is greater than the maximum theoretical dE
+ //(dE > GetMaxDeltaE(Z,A)) then we return the incident energy corresponding to the maximum,
+ // GetEIncOfMaxDeltaE(Z,A)
+
+ if(Z<1) return 0.;
+
+ Double_t DE = (delta_e > 0 ? delta_e : GetEnergyLoss());
+
+ TF1* dE = GetELossFunction(Z, A);
+ Double_t e1,e2;
+ dE->GetRange(e1,e2);
+ switch(type){
+ case KVIonRangeTable::kEmin:
+ e2=GetEIncOfMaxDeltaE(Z,A);
+ break;
+ case KVIonRangeTable::kEmax:
+ e1=GetEIncOfMaxDeltaE(Z,A);
+ break;
+ }
+ return dE->GetX(DE,e1,e2);
+}
+
+Double_t KVDetector::GetDeltaEFromERes(Int_t Z, Int_t A, Double_t Eres)
+{
+ // Overrides KVMaterial::GetDeltaEFromERes
+ //
+ // Calculate energy loss in active layer of detector for nucleus (Z,A)
+ // having a residual kinetic energy Eres (MeV)
+
+ if(Z<1 || Eres<KVDETECTOR_MINIMUM_E) return 0.;
+ Double_t Einc = GetIncidentEnergyFromERes(Z,A,Eres);
+ if(Einc < KVDETECTOR_MINIMUM_E) return 0.;
+ return GetELossFunction(Z,A)->Eval(Einc);
+}
+
+Double_t KVDetector::GetIncidentEnergyFromERes(Int_t Z, Int_t A, Double_t Eres)
+{
+ // Overrides KVMaterial::GetIncidentEnergyFromERes
+ //
+ // Calculate incident energy of nucleus from residual energy
+
+ if(Z<1 || Eres<KVDETECTOR_MINIMUM_E) return 0.;
+ return GetEResFunction(Z,A)->GetX(Eres);
+}
+
+Double_t KVDetector::GetSmallestEmaxValid(Int_t Z, Int_t A)
+{
+ // Returns the smallest maximum energy for which range tables are valid
+ // for all absorbers in the detector, and given ion (Z,A)
+
+ Double_t maxmin=-1.;
+ TIter next(fAbsorbers);
+ KVMaterial *abs;
+ while( (abs = (KVMaterial*)next()) ){
+ if(maxmin<0.) maxmin = abs->GetEmaxValid(Z,A);
+ else {
+ if(abs->GetEmaxValid(Z,A) < maxmin) maxmin = abs->GetEmaxValid(Z,A);
+ }
+ }
+ return maxmin;
+}
+
+void KVDetector::ReadDefinitionFromFile(const Char_t* envrc)
+{
+ // Create detector from text file in 'TEnv' format.
+ //
+ // Example:
+ // ========
+ //
+ // Layer: Gold
+ // Gold.Material: Au
+ // Gold.AreaDensity: 200.*KVUnits::ug
+ // +Layer: Gas1
+ // Gas1.Material: C3F8
+ // Gas1.Thickness: 5.*KVUnits::cm
+ // Gas1.Pressure: 50.*KVUnits::mbar
+ // Gas1.Active: yes
+ // +Layer: Si1
+ // Si1.Material: Si
+ // Si1.Thickness: 300.*KVUnits::um
+
+ TEnv fEnvFile(envrc);
+
+ KVString layers(fEnvFile.GetValue("Layer", ""));
+ layers.Begin(" ");
+ while( !layers.End() ){
+ KVString layer = layers.Next();
+ KVString mat = fEnvFile.GetValue(Form("%s.Material",layer.Data()),"");
+ KVString tS = fEnvFile.GetValue(Form("%s.Thickness",layer.Data()),"");
+ KVString pS = fEnvFile.GetValue(Form("%s.Pressure",layer.Data()),"");
+ KVString dS = fEnvFile.GetValue(Form("%s.AreaDensity",layer.Data()),"");
+ Double_t thick, dens, press;
+ thick=dens=press=0.;KVMaterial* M = 0;
+ if(pS != "" && tS != ""){
+ press = (Double_t)gROOT->ProcessLineFast( Form("%s*1.e+12", pS.Data()) );
+ press/=1.e+12;
+ thick = (Double_t)gROOT->ProcessLineFast( Form("%s*1.e+12", tS.Data()) );
+ thick/=1.e+12;
+ M = new KVMaterial(mat.Data(), thick, press);
+ }
+ else if(tS != ""){
+ thick = (Double_t)gROOT->ProcessLineFast( Form("%s*1.e+12", tS.Data()) );
+ thick/=1.e+12;
+ M = new KVMaterial(mat.Data(), thick);
+ }
+ else if(dS != ""){
+ dens = (Double_t)gROOT->ProcessLineFast( Form("%s*1.e+12", dS.Data()) );
+ dens/=1.e+12;
+ M = new KVMaterial(dens, mat.Data());
+ }
+ if(M){
+ AddAbsorber(M);
+ if(fEnvFile.GetValue(Form("%s.Active",layer.Data()),kFALSE)) SetActiveLayer(M);
+ }
+ }
+}
=== modified file 'KVMultiDet/detectors/KVDetector.h'
--- KVMultiDet/detectors/KVDetector.h 2010-12-03 09:16:52 +0000
+++ KVMultiDet/detectors/KVDetector.h 2011-03-14 15:58:00 +0000
@@ -27,12 +27,6 @@
#define KVD_NORECPRC_CNXN 0
#endif
-
-#define KVDETECTOR_UNKNOWN_MATERIAL "Constructor called with unknown material"
-#define KVDETECTOR_UNKNOWN_DETECTOR "Constructor called with unknown detector prototype"
-#define KVDETECTOR_NOT_IN_TELESCOPE "Detector not associated to any telescope: position undefined"
-#define KVDETECTOR_ADD_TO_UNKNOWN_TELESCOPE "Pointer to telescope is null."
-
//for energies less than this (MeV) particles are considered to be stopped
#define KVDETECTOR_MINIMUM_E 0.001
@@ -51,9 +45,6 @@
class TGeoVolume;
class TTree;
-Double_t ELossActive(Double_t * x, Double_t * par);
-Double_t EResDet(Double_t * x, Double_t * par);
-
class KVDetector:public KVMaterial {
friend class KVDetectorBrowser;
@@ -78,10 +69,6 @@
Int_t fUnidentP; //! temporary counters, determine state of identified/unidentified particle flags
protected:
-Int_t npar_loss;//!number of params for eloss function
-Int_t npar_res;//!number of params for eres function
-Double_t *par_loss;//!array of params for eloss function
-Double_t *par_res;//!array of params for eres function
TString fFName; //!dynamically generated full name of detector
KVList *fModules; //references to connected electronic modules (not implemented yet)
@@ -98,6 +85,14 @@
Binary8_t fFiredMask;//bitmask used by Fired to determine which parameters to take into account
+ Double_t ELossActive(Double_t * x, Double_t * par);
+ Double_t EResDet(Double_t * x, Double_t * par);
+
+ TF1* fELossF; //! parametric function dE in active layer vs. incident energy
+ TF1* fEResF; //! parametric function Eres residual energy after all layers of detector
+
+ Double_t fEResforEinc;//! used by GetIncidentEnergy & GetCorrectedEnergy
+
public:
KVDetector();
KVDetector(const Char_t * type, const Float_t thick = 0.0);
@@ -129,14 +124,12 @@
Double_t GetTotalThicknessInCM()
{
- // Calculate and return the total thickness of ALL absorbers making up the detector,
- // not just the active layer (as for GetThickness()).
- //
- // As different materials can have different default units for thickness, we convert
- // everything into centimetres and return the total thickness with these units.
+ // Calculate and return the total thickness in centimetres of ALL absorbers making up the detector,
+ // not just the active layer (value returned by GetThickness()).
+
fTotThickness=0;
TIter next(fAbsorbers); KVMaterial* mat;
- while( (mat = (KVMaterial*)next()) ) fTotThickness += mat->GetThicknessInCM();
+ while( (mat = (KVMaterial*)next()) ) fTotThickness += mat->GetThickness();
return fTotThickness;
};
@@ -171,7 +164,7 @@
};
virtual Double_t GetCorrectedEnergy(UInt_t z, UInt_t a, Double_t e =
-1., Bool_t transmission=kTRUE);
- virtual UInt_t FindZmin(Double_t ELOSS = -1., Char_t mass_formula = -1);
+ virtual Int_t FindZmin(Double_t ELOSS = -1., Char_t mass_formula = -1);
void AddACQParam(KVACQParam*);
virtual KVACQParam *GetACQParam(const Char_t*/*name*/);
@@ -264,9 +257,6 @@
return TestBit(kIsBeingDeleted);
};
- virtual void SetEResParams(Int_t Z, Int_t A);
- virtual void SetELossParams(Int_t Z, Int_t A);
-
static KVDetector *MakeDetector(const Char_t * name, Float_t thick);
const TVector3& GetNormal();
@@ -287,6 +277,31 @@
// detector implementations: this version just returns -1.
return -1;
};
+ virtual Double_t GetEIncOfMaxDeltaE(Int_t Z, Int_t A);
+ virtual Double_t GetDeltaE(Int_t Z, Int_t A, Double_t Einc);
+ virtual Double_t GetTotalDeltaE(Int_t Z, Int_t A, Double_t Einc);
+ virtual Double_t GetERes(Int_t Z, Int_t A, Double_t Einc);
+ virtual Double_t GetIncidentEnergy(Int_t Z, Int_t A, Double_t delta_e =
+ -1.0, enum KVIonRangeTable::SolType type = KVIonRangeTable::kEmax);
+ /*virtual Double_t GetEResFromDeltaE(...) - DON'T IMPLEMENT, CALLS GETINCIDENTENERGY*/
+ virtual Double_t GetDeltaEFromERes(Int_t Z, Int_t A, Double_t Eres);
+ virtual Double_t GetIncidentEnergyFromERes(Int_t Z, Int_t A, Double_t Eres);
+
+ virtual TF1* GetEResFunction(Int_t Z, Int_t A);
+ virtual TF1* GetELossFunction(Int_t Z, Int_t A);
+
+ virtual Double_t GetSmallestEmaxValid(Int_t Z, Int_t A);
+
+ virtual void SetEResAfterDetector(Double_t e)
+ {
+ fEResforEinc=e;
+ };
+ virtual Double_t GetEResAfterDetector() const
+ {
+ return fEResforEinc;
+ };
+
+ virtual void ReadDefinitionFromFile(const Char_t*);
ClassDef(KVDetector, 8) //Base class for the description of detectors in multidetector arrays
};
=== modified file 'KVMultiDet/detectors/KVMaterial.cpp'
--- KVMultiDet/detectors/KVMaterial.cpp 2009-11-26 16:48:10 +0000
+++ KVMultiDet/detectors/KVMaterial.cpp 2011-03-14 15:58:00 +0000
@@ -1,5 +1,4 @@
/***************************************************************************
-$Id: KVMaterial.cpp,v 1.55 2009/04/01 09:28:02 franklan Exp $
kvmaterial.cpp - description
-------------------
begin : Thu May 16 2002
@@ -33,65 +32,31 @@
///////////////////////////////////////////////////////////////////////////////
// KVMaterial
//
-//The KVMaterial class encapsulates the routines of R. Dayras et al for calculating
-//the energy lost by charged particles (nuclei) in different absorbers,
-//formerly known as the "vedaloss" routines.
-//
-//To create an absorber, use the ctor specifying a type of material and optionally
-//the thickness of the absorber.
-// KVMaterial mat1("Si",150.0); //150 microns of silicon
-// KVMaterial mat2("Au",0.8); //800 micrograms/cm**2 of gold
-//The available materials are defined in the $KVROOT/KVFiles/kvloss.data file. In order to get
-//a list of the currently defined materials, use KVMaterial::GetListOfMaterials().
-//
-//Depending on the type of material, the default units of "thickness" are different:
-// fType = Si, Myl - thickness in micrometres
-// fType = NE102, CsI - thickness in centimetres
-// fType = C3F8, CF4 - pressure in mbar
-// fType = Au, Sn, U, Nb, Ni, C, Ta - thickness in mg/cm2
-//
-//By default, the naturally-occuring atomic mass is used for the elements. To use
-//some rarer isotope, use the SetMass() method with the appropriate mass (not
-//valid for compounds).
-//
+// Description of materials used for detectors and targets.
///////////////////////////////////////////////////////////////////////////////
-TString KVMaterial::kDataFilePath = "";
-Char_t KVMaterial::kUnits[][10] = {
- "mbar",
- "microns",
- "mg/cm2",
- "cm",
- "Torr"
-};
+KVIonRangeTable* KVMaterial::fIonRangeTable = 0x0;
//___________________________________________________________________________________
void KVMaterial::init()
{
- //Default initialisations.
- //No properties are set for the material.
- fAmasr = 0.;
- fAmat = 0.;
- fDens = 0.;
- fGasThick = 0.;
- fMoleWt = 0.;
- fTemp = 19.;
- fUnits = 10;
- fZmat = 0.;
- if (kDataFilePath == "") {
- //get full path to range tables file from env vars (i.e. .kvrootrc)
- kDataFilePath = gEnv->GetValue("KVMaterial.RangeTables", "");
- gSystem->ExpandPathName( kDataFilePath );
- }
- for (int i = 0; i < 100; i++) {
- fEmin[i] = 0.0;
- fEmax[i] = 500.0;
- }
+ // Default initialisations.
+ // No properties are set for the material (except standard temperature (19°C) and pressure (1 atm))
+ // Default range table is generated if not already done.
+ // By default it is the VEDALOSS table implemented in KVedaLoss.
+ // You can change this by changing the value of environment variable KVMaterial.IonRangeTable.
+
fELoss = 0;
SetName("");
SetTitle("");
- ELoss = ERes = 0;
fNormToMat.SetXYZ(0,0,1);
+ fAmasr = 0;
+ fPressure = 1. * KVUnits::atm;
+ fTemp = 19.0;
+ // create default range table singleton if not already done
+ if(!fIonRangeTable) {
+ fIonRangeTable = KVIonRangeTable::GetRangeTable( gEnv->GetValue("KVMaterial.IonRangeTable", "VEDALOSS") );
+ }
}
//
@@ -102,11 +67,37 @@
}
//__________________________________________________________________________________
-KVMaterial::KVMaterial(const Char_t * type, const Float_t thick)
-{
- //Initialise absorber with given type of material and thickness (default = 0.0)
- init();
- SetMaterial(type);
+KVMaterial::KVMaterial(const Char_t * type, const Double_t thick)
+{
+ // Create material with given type and linear thickness in cm.
+
+ init();
+ SetMaterial(type);
+ SetThickness(thick);
+}
+
+KVMaterial::KVMaterial(Double_t area_density, const Char_t * type)
+{
+ // Create material with given area density in g/cm**2 and given type
+
+ init();
+ SetMaterial(type);
+ SetAreaDensity(area_density);
+}
+
+KVMaterial::KVMaterial(const Char_t * gas, const Double_t thick, const Double_t pressure, const Double_t temperature)
+{
+ // Create gaseous material with given type, linear thickness in cm, pressure in Torr,
+ // and temperature in degrees C (default value 19°C).
+ //
+ // Examples
+ // 15 cm of CF4 gas at 1 atm and 19°C : KVMaterial("CF4", 15., 1.*KVUnits::atm)
+ // 50 mm of C3F8 at 30 mbar and 25°C : KVMaterial("C3F8", 50.*KVUnits::mm, 30.*KVUnits::mbar, 25.)
+
+ init();
+ SetMaterial(gas);
+ fPressure = pressure;
+ fTemp = temperature;
SetThickness(thick);
}
@@ -123,77 +114,14 @@
}
//___________________________________________________________________________________
-KVList *KVMaterial::GetListOfMaterials()
-{
- //Creates and fills a list with all the materials currently described in the
- //$KVROOT/KVFiles/kvloss.data parameters file.
- //To avoid memory leakage, delete the list after use !!
-
- KVList *tmp = new KVList();
-
- Char_t name[25], gtype[25];
- Float_t Amat = 0.;
- Float_t Dens = 0.;
- Float_t GasThick = 0.;
- Float_t MoleWt = 0.;
- Float_t Temp = 19.;
- UInt_t Units = 10;
- Float_t Zmat = 0.;
-
- if (kDataFilePath == "") {
- //get full path to range tables file from env vars (i.e. .kvrootrc)
- kDataFilePath = gEnv->GetValue("KVMaterial.RangeTables", "");
- gSystem->ExpandPathName( kDataFilePath );
- }
-
- FILE *fp;
- if (!(fp = fopen(kDataFilePath.Data(), "r"))) {
- printf("KVMaterial::GetMaterialsList : ");
- printf(KVMATERIAL_INIT_FILE_NOT_FOUND, kDataFilePath.Data());
- return 0;
- } else {
- char line[132];
- while (fgets(line, 132, fp)) { // read lines from file
-
- switch (line[0]) {
-
- case '/': // ignore comment lines
- break;
-
- case '+': // header lines
-
- if (sscanf(line, "+ %s %s %u %f %f %f %f %f %f",
- gtype, name, &Units, &Dens, &Zmat, &Amat,
- &GasThick, &MoleWt, &Temp)
- != 9) {
- cout << "KVMaterial::GetMaterialsList> : ";
- printf(KVMATERIAL_INIT_FILE_HEADER_PROBLEM, kDataFilePath.Data());
- cout << endl;
- fclose(fp);
- return 0;
- }
-//found a new material
- KVMaterial *tmp_mat = new KVMaterial(gtype);
- tmp->Add(tmp_mat);
- break;
- }
- }
- fclose(fp);
- }
- return tmp;
-}
-
-//___________________________________________________________________________________
void KVMaterial::SetMaterial(const Char_t * mat_type)
{
- //Intialise material of a given type.
- //The known types are defined by the contents of file $KVROOT/kvloss.data.
- //For materials which are elements of the periodic table you can specify
- //the isotope such as "64Ni", "13C", "natSn", etc. etc.
+ // Intialise material of a given type.
+ // The material must exist in the currently used range tables (fIonRangeTable).
+ // For materials which are elements of the periodic table you can specify
+ // the isotope such as "64Ni", "13C", "natSn", etc. etc.
- Char_t name[25], gtype[25];
init();
- FILE *fp;
//are we dealing with an isotope ?
Char_t type[10];
Int_t iso_mass = 0;
@@ -202,116 +130,29 @@
strcpy(type, mat_type);
}
}
-
- if (!(fp = fopen(kDataFilePath.Data(), "r"))) {
- Error("SetMaterial", KVMATERIAL_INIT_FILE_NOT_FOUND,
- kDataFilePath.Data());
- } else {
- char line[132];
- while (fgets(line, 132, fp)) { // read lines from file
-
- switch (line[0]) {
-
- case '/': // ignore comment lines
- break;
-
- case '+': // header lines
-
- if (sscanf(line, "+ %s %s %u %lf %lf %lf %lf %lf %lf",
- gtype, name, &fUnits, &fDens, &fZmat, &fAmat,
- &fGasThick, &fMoleWt, &fTemp)
- != 9) {
- Error("SetMaterial", KVMATERIAL_INIT_FILE_HEADER_PROBLEM,
- kDataFilePath.Data());
- goto bad_exit;
- }
- SetName(name);
- SetType(gtype);
- if (!strcmp(GetType(), type)) { // is this what you're looking for ?
- //look for energy limits to calculation validity
- if (!fgets(line, 132, fp)) {
- Warning("SetMaterial", KVMATERIAL_INIT_FILE_DATA_PROBLEM,
- kDataFilePath.Data());
- goto bad_exit;
- } else {
- while (line[0] == 'Z') {
- Int_t z1, z2;
- Float_t e1, e2;
- sscanf(line, "Z = %d,%d %f < E/A < %f MeV", &z1,
- &z2, &e1, &e2);
- for (int i = z1; i <= z2; i++) {
- fEmin[i - 1] = e1;
- fEmax[i - 1] = e2;
- }
- fgets(line, 132, fp);
- }
- }
-
- for (register int count = 0; count < 100; count++) {
-
- if (sscanf(line, "%lf %lf %lf %lf %lf %lf %lf %lf",
- &fCoeff[count][0], &fCoeff[count][1],
- &fCoeff[count][2], &fCoeff[count][3],
- &fCoeff[count][4], &fCoeff[count][5],
- &fCoeff[count][6], &fCoeff[count][7])
- != 8) {
- Error("SetMaterial",
- KVMATERIAL_INIT_FILE_COEFF_PROBLEM, kDataFilePath.Data(),
- type);
- goto bad_exit;
- }
- if (!fgets(line, 132, fp)) {
- Warning("SetMaterial",
- KVMATERIAL_INIT_FILE_DATA_PROBLEM,
- kDataFilePath.Data());
- goto bad_exit;
- } else {
- if (sscanf(line, "%lf %lf %lf %lf %lf %lf",
- &fCoeff[count][8], &fCoeff[count][9],
- &fCoeff[count][10], &fCoeff[count][11],
- &fCoeff[count][12], &fCoeff[count][13])
- != 6) {
- Error("SetMaterial",
- KVMATERIAL_INIT_FILE_COEFF_PROBLEM,
- kDataFilePath.Data(), type);
- goto bad_exit;
- }
- }
- fgets(line, 132, fp);
- }
- goto normal_exit; // i apologise for the use of "g**o"
- }
- break;
- }
- }
-// only come here if nothing has been found
- Warning("SetMaterial", KVMATERIAL_TYPE_NOT_FOUND, type, kDataFilePath.Data());
- init(); //no properties are set
- return;
- bad_exit:
- Fatal("SetMaterial", "Programme execution impossible");
- normal_exit:
- fclose(fp);
- //set isotope in case call was made with "64Ni" or "48Ca" etc.
- if (iso_mass)
- SetMass(iso_mass);
- }
+ if(iso_mass) SetMass( iso_mass );
+ SetType(type);
+ if(!fIonRangeTable->IsMaterialKnown(type))
+ Warning("SetMaterial",
+ "Called for material %s which is unknown in current range table %s. Energy loss & range calculations impossible.",
+ type, fIonRangeTable->GetName());
+ else {
+ SetType( fIonRangeTable->GetMaterialName(type) );
+ SetName(mat_type);
+ }
}
//___________________________________________________________________________________
KVMaterial::~KVMaterial()
{
//Destructor
- //The TF1 objects "ELoss" and "ERes" are not to be deleted as they can be
- //used by other KVMaterial objects (they are retrieved from gROOT->GetListOfFunctions())
- //and cannot cause a memory leak as they are only ever created if not found in
- //the gROOT list of functions.
}
-void KVMaterial::SetMass(Float_t a)
+void KVMaterial::SetMass(Double_t a)
{
-//Set the atomic mass of the material - use if you want to change the default naturally
-//occuring mass for some rarer isotope.
+ //Set the atomic mass of the material - use if you want to change the default naturally
+ //occuring mass for some rarer isotope.
+
if (GetActiveLayer()) {
GetActiveLayer()->SetMass(a);
return;
@@ -323,9 +164,10 @@
Double_t KVMaterial::GetMass() const
{
//Returns atomic mass of material. Will be isotopic mass if set.
+
if (GetActiveLayer())
return GetActiveLayer()->GetMass();
- return (fAmasr ? fAmasr : fAmat);
+ return (fAmasr ? fAmasr : fIonRangeTable->GetAtomicMass(GetType()));
}
//___________________________________________________________________________________
@@ -353,81 +195,154 @@
}
//___________________________________________________________________________________
+
+Bool_t KVMaterial::IsGas() const
+{
+ // Returns kTRUE for gaseous materials/detectors.
+
+ if (GetActiveLayer())
+ return GetActiveLayer()->IsGas();
+ return fIonRangeTable->IsMaterialGas(GetType());
+}
+
+//___________________________________________________________________________________
+
Double_t KVMaterial::GetZ() const
{
//Returns atomic number of material.
if (GetActiveLayer())
return GetActiveLayer()->GetZ();
- return fZmat;
+ return fIonRangeTable->GetZ(GetType());
}
//___________________________________________________________________________________
+
Double_t KVMaterial::GetDensity() const
{
- //Returns density of material.
- //For solids, the units are g/cm**3.
- //For a gas, density is calculated from molecular weight, pressure
- //and temperature according to ideal gas law
+ // Returns density of material in g/cm**3.
+ // For a gas, density is calculated from current pressure & temperature according to ideal gas law
if (GetActiveLayer())
return GetActiveLayer()->GetDensity();
- if (fUnits == kTORR)
- return (fMoleWt * fThick) / ((fTemp + ZERO_KELVIN) * RTT);
- if (fUnits == kMBAR)
- return (fMoleWt * fThick * 0.75) / ((fTemp + ZERO_KELVIN) * RTT);
- return fDens;
-}
-
-void KVMaterial::SetPressure(Float_t t)
-{
- //Set the pressure of a gas.
- //If the material is not a gas, this has no effect.
- //The units depend on fUnits
-
- if (GetActiveLayer()) {
- GetActiveLayer()->SetPressure(t);
- return;
- }
-
- if (fUnits == kMBAR || fUnits == kTORR)
- fThick = t;
-}
-
-
-//___________________________________________________________________________________
-
-Float_t KVMaterial::GetPressure() const
-{
- //Returns the pressure of a gas.
- //If the material is not a gas, value is zero.
- //The units depend on fUnits
-
+ fIonRangeTable->SetTemperatureAndPressure(GetType(), fTemp, fPressure);
+ return fIonRangeTable->GetDensity(GetType());
+}
+
+//___________________________________________________________________________________
+
+void KVMaterial::SetThickness(Double_t t)
+{
+ // Set the linear thickness of the material in cm or use one of the
+ // Units constants:
+ // SetThickness( 30.*KVUnits::um ); set thickness to 30 microns
+
+ if (GetActiveLayer()) {
+ GetActiveLayer()->SetThickness(t);
+ return;
+ }
+ // recalculate area density
+ fThick = t * GetDensity();
+}
+
+//___________________________________________________________________________________
+
+Double_t KVMaterial::GetThickness() const
+{
+ // Returns the linear thickness of the material in cm.
+ // Use Units to change units:
+ // mat.GetThickness()/KVUnits::um ; in microns
+
+ if (GetActiveLayer())
+ return GetActiveLayer()->GetThickness();
+ return fThick/GetDensity();
+}
+
+//___________________________________________________________________________________
+
+void KVMaterial::SetAreaDensity(Double_t dens /* g/cm**2 */)
+{
+ // Set area density in g/cm**2.
+ //
+ // For solids, area density can only changed by changing linear dimension
+ // (fixed density).
+ // For gases, the density depends on temperature and pressure. This method
+ // leaves temperature and pressure unchanged, therefore for gases also this
+ // method will effectively modify the linear dimension of the gas cell.
+
+ if(GetActiveLayer())
+ GetActiveLayer()->SetAreaDensity(dens);
+ fThick = dens;
+}
+
+//___________________________________________________________________________________
+
+Double_t KVMaterial::GetAreaDensity() const
+{
+ // Return area density of material in g/cm**2
+
+ if(GetActiveLayer()) return GetActiveLayer()->GetAreaDensity();
+ return fThick;
+}
+
+//___________________________________________________________________________________
+
+void KVMaterial::SetPressure(Double_t p)
+{
+ // Set the pressure of a gaseous material (in torr)
+ // As this changes the density of the gas, it also changes
+ // the area density of the absorber (for fixed linear dimension)
+
+ if(!IsGas()) return;
+ if (GetActiveLayer()) {
+ GetActiveLayer()->SetPressure(p);
+ return;
+ }
+ // get current linear dimension of absorber
+ Double_t e = GetThickness();
+ fPressure = p;
+ // change area density to keep linear dimension constant
+ SetThickness(e);
+}
+
+
+//___________________________________________________________________________________
+
+Double_t KVMaterial::GetPressure() const
+{
+ // Returns the pressure of a gas (in torr).
+ // If the material is not a gas, value is zero.
+
+ if(!IsGas()) return 0.0;
if (GetActiveLayer())
return GetActiveLayer()->GetPressure();
- if (fUnits == kMBAR || fUnits == kTORR)
- return fThick;
- return 0.0;
+ return fPressure;
}
//___________________________________________________________________________________
-void KVMaterial::SetTemperature(Float_t t)
+void KVMaterial::SetTemperature(Double_t t)
{
- //Set temperature of material.
- //The units are: degrees celsius
+ // Set temperature of material.
+ // The units are: degrees celsius
+ // As this changes the density of the gas, it also changes
+ // the area density of the absorber (for fixed linear dimension)
+ if(!IsGas()) return;
if (GetActiveLayer()) {
GetActiveLayer()->SetTemperature(t);
return;
}
-
+ // get current linear dimension of absorber
+ Double_t e = GetThickness();
fTemp = t;
+ // change area density to keep linear dimension constant
+ SetThickness(e);
}
//___________________________________________________________________________________
-Float_t KVMaterial::GetTemperature() const
+Double_t KVMaterial::GetTemperature() const
{
//Returns temperature of material.
//The units are: degrees celsius
@@ -439,105 +354,11 @@
//___________________________________________________________________________________
-void KVMaterial::SetUnits(UInt_t u)
-{
- //Set units for 'thickness' of material.
- //WARNING: changing the default units will alter how the 'thickness'
- //is used in energy loss calculations!
- //WARNING2: we do not perform a conversion of any existing value
- //for the 'thickness' of the material, the existing value will be used
- //with the new units
-
- if (GetActiveLayer()) {
- GetActiveLayer()->SetUnits(u);
- return;
- }
-
- fUnits = u;
-}
-
-
-//___________________________________________________________________________________
-
-UInt_t KVMaterial::GetUnits() const
-{
- //Returns units for 'thickness' of material.
- //Can be one of KVMaterial::kMBAR, KVMaterial::kMGCM2, KVMaterial::kTORR,
- //KVMaterial::kMICRON, KVMaterial::kCM.
-
- if (GetActiveLayer())
- return GetActiveLayer()->GetUnits();
- return fUnits;
-}
-
-//___________________________________________________________________________________
-
-void KVMaterial::SetThickness(Float_t t)
-{
- //Set the thickness of the material.
- //
- //The meaning and units of "thickness" depend on the material.
- //For a gas it is the depth of the gas cell in millimetres.
- //For silicon and mylar it is the thickness in micrometres.
- //For CsI and plastic scintillators it is the length/thickness in centimetres.
- //For metallic elements (such as for targets) it is the areal density
- //in milligrams per cm^2.
-
- if (GetActiveLayer()) {
- GetActiveLayer()->SetThickness(t);
- return;
- }
-
- if (fUnits == kMBAR || fUnits == kTORR)
- fGasThick = t; //thickness of gas cell in cm
- else
- fThick = t;
-}
-
-//___________________________________________________________________________________
-
-Float_t KVMaterial::GetThickness() const
-{
- //Returns the "thickness" of the material - units depend on the material type
- //For a gas it is the depth of the gas cell in millimetres.
- //For silicon and mylar it is the thickness in micrometres.
- //For CsI and plastic scintillators it is the length/thickness in centimetres.
- //For metallic elements (such as for targets) it is the areal density
- //in milligrams per cm^2.
-
- if (GetActiveLayer())
- return GetActiveLayer()->GetThickness();
- if (fUnits == kMBAR || fUnits == kTORR)
- return fGasThick;
- return fThick;
-}
-
-//___________________________________________________________________________________
-
-Double_t KVMaterial::GetThicknessInCM() const
-{
- //Returns the thickness of the material in CENTIMETRES
-
- if (GetActiveLayer())
- return GetActiveLayer()->GetThicknessInCM();
- Double_t t = (Double_t)fGasThick;
- if (fUnits == kMBAR || fUnits == kTORR)
- return (t/10.);
- t = (Double_t)fThick;
- if(fUnits == kMICRONS)
- return (t/10000.);
- else if(fUnits == kMGCM2)
- return (t/fDens);
- return fThick;
-}
-
-//___________________________________________________________________________________
-
Double_t KVMaterial::GetEffectiveThickness(TVector3 & norm,
TVector3 & direction)
{
- //Calculate effective thickness of absorber as 'seen' in 'direction', taking into
- //account the arbitrary orientation of the 'norm' normal to the material's surface
+ // Calculate effective linear thickness of absorber (in cm) as 'seen' in 'direction', taking into
+ // account the arbitrary orientation of the 'norm' normal to the material's surface
TVector3 n = norm.Unit();
TVector3 d = direction.Unit();
@@ -547,248 +368,57 @@
}
//___________________________________________________________________________________
+
+Double_t KVMaterial::GetEffectiveAreaDensity(TVector3 & norm,
+ TVector3 & direction)
+{
+ // Calculate effective area density of absorber (in g/cm**2) as 'seen' in 'direction', taking into
+ // account the arbitrary orientation of the 'norm' normal to the material's surface
+
+ TVector3 n = norm.Unit();
+ TVector3 d = direction.Unit();
+ //absolute value of scalar product, in case direction is opposite to normal
+ Double_t prod = TMath::Abs(n * d);
+ return GetAreaDensity() / TMath::Max(prod, 1.e-03);
+}
+
+//___________________________________________________________________________________
void KVMaterial::Print(Option_t * option) const
{
//Show information on this material
cout << "KVMaterial: " << GetName() << " (" << GetType() << ")" << endl;
- if (fUnits == kMBAR || fUnits == kTORR)
- cout << " Thickness " << KVMaterial::
- GetThickness() << " mm, pressure " << GetPressure() << " " <<
- GetThicknessUnits() << endl;
- else
- cout << " Thickness " << KVMaterial::
- GetThickness() << " " << GetThicknessUnits() << endl;
+ if (fIonRangeTable->IsMaterialGas(GetType()))
+ cout << " Pressure " << GetPressure() << " torr" << endl;
+ cout << " Thickness " << KVMaterial::GetThickness() << " cm" << endl;
+ cout << " Area density " << KVMaterial::GetAreaDensity() << " g/cm**2" << endl;
cout << "-----------------------------------------------" << endl;
cout << " Z = " << GetZ() << " atomic mass = " << GetMass() << endl;
- cout << " Density = " << GetDensity() << endl;
- cout << " Validity of calculation:" << endl;
- cout << " Z = 1,2 -- " << GetEminVedaloss(1) << " AMeV < E < " <<
- GetEmaxVedaloss(1) << " AMeV" << endl;
- cout << " Z >= 3 -- " << GetEminVedaloss(3) << " AMeV < E < " <<
- GetEmaxVedaloss(3) << " AMeV" << endl;
+ cout << " Density = " << GetDensity() << " g/cm**3" << endl;
cout << "-----------------------------------------------" << endl;
}
-//_______________________________________________________________________________________
-Double_t ELossSaclay(Double_t * x, Double_t * par)
-{
- //Calculates energy loss (DE) as a function of incident energy x[0]
- //Parameters:
- // par[0] - par[14] : Saclay range tables parameters
- // par[0] = Z of nucleus
- // par[15] = A of nucleus
- // par[16] = minimum incident energy in MeV for valid calculation
- // par[17] = maximum incident energy in MeV for valid calculation
- // par[18] = ratio of isotopic to natural mass of material, log(fAmasr/fAmat)
-
- //energy loss = incident energy - residual energy
- return x[0] - EResSaclay(x, par);
-}
-
-//_______________________________________________________________________________________
-Double_t EResSaclay(Double_t * x, Double_t * par)
-{
- //Calculates residual energy after passage through material
- //as a function of incident energy x[0]
- //Parameters:
- // par[0] - par[14] : Saclay range tables parameters
- // par[0] = Z of nucleus
- // par[15] = A of nucleus
- // par[16] = minimum incident energy in MeV for valid calculation
- // par[17] = maximum incident energy in MeV for valid calculation
- // par[18] = ratio of isotopic to natural mass of material, log(fAmasr/fAmat)
-
- //check incident energy is valid
- //incident energy less than minimum valid energy - residual energy=0 (particle stops)
- //if (x[0] < par[16])
- // return 0.0;
- //incident energy greater than maximum valid energy - residual E = incident E
- //if (x[0] > par[17])
- // return x[0];
-
- // set up polynomial
- Double_t x1 = TMath::Log(0.1);
- Double_t x2 = TMath::Log(0.2);
- Double_t ran = 0.0;
- for (register int j = 2; j < 7; j++)
- ran += par[j + 1] * TMath::Power(x2, (Double_t) (j - 1));
- ran += par[2];
- Double_t y2 = ran;
- ran = 0.0;
- for (register int jj = 2; jj < 7; jj++)
- ran += par[jj + 1] * TMath::Power(x1, (Double_t) (jj - 1));
- ran += par[2];
- Double_t y1 = ran;
- Double_t adm = (y2 - y1) / (x2 - x1);
- Double_t adn = (y1 - adm * x1);
- Double_t arm = y1;
- // calculate energy loss
- Double_t eps = x[0] / par[15]; //energy in MeV/nucleon
- Double_t dleps = TMath::Log(eps);
- Double_t riso = TMath::Log(par[15] / par[1]) + par[18];
-
- if (eps < 0.1)
- ran = adm * dleps + adn;
- else {
- ran = 0.0;
- for (register int j = 2; j < 7; j++)
- ran += par[j + 1] * TMath::Power(dleps, (Double_t) (j - 1));
- ran += par[2];
- }
- ran += riso;
-
- Double_t range = TMath::Exp(ran);
- range -= par[14];
- if (range <= 0.005) { // particle stopped in material, gives up all energy
- return 0.0; //residual energy = 0
- }
-
- Double_t ranx = TMath::Log(range);
- Double_t ranx1 = ranx - riso;
- Double_t depsx;
- if (ranx1 < arm)
- depsx = (ranx1 - adn) / adm;
- else {
- depsx = 0.0;
- for (register int j = 2; j < 7; j++)
- depsx += par[j + 7] * TMath::Power(ranx1, (Double_t) (j - 1));
- depsx += par[8];
- }
-
- Double_t eps1 = depsx + TMath::Log(1 - PERC);
- Double_t eps2 = depsx + TMath::Log(1 + PERC);
- Double_t rap = TMath::Log((1 + PERC) / (1 - PERC));
-
- Double_t rn1 = 0.0;
- if (TMath::Exp(eps1) < 0.1)
- rn1 = adm * eps1 + adn;
- else {
- for (register int j = 1; j < 7; j++)
- rn1 += par[j + 1] * TMath::Power(eps1, (Double_t) (j - 1));
- }
- Double_t rn2 = 0.0;
- if (TMath::Exp(eps2) < 0.1)
- rn2 = adm * eps2 + adn;
- else {
- for (register int j = 1; j < 7; j++)
- rn2 += par[j + 1] * TMath::Power(eps2, (Double_t) (j - 1));
- }
-
- Double_t epres = eps1 + (rap / (rn2 - rn1)) * (ranx1 - rn1);
- epres = TMath::Exp(epres);
- Double_t eres = par[15] * epres;
- //make sure residual energy is <= incident energy
- //following original code, the minimum energy loss is 0.005 MeV
-
- Double_t eloss = x[0] - eres;
- if (eloss < 0.005 && eres > eloss)
- eres = x[0];
- return eres;
-}
-
-//_______________________________________________________________________________________//
-
-void KVMaterial::GetELossParams(Int_t Z, Int_t A, Double_t * par)
-{
- // get coefficients associated with this nucleus Z
- if(Z<1 || Z>ZMAX_VEDALOSS) return;
-
- for (int i = 0; i < 14; i++) {
-
- par[i] = fCoeff[(Z - 1)][i];
-
- }
- fAmasr = (fAmasr ? fAmasr : fAmat);
- // calculate "thickness"
- //for isotopic materials, the density used in the conversion
- //from microns or cm to mg/cm2 is modified according to the
- //chosen isotopic mass, fAmasr
- switch (fUnits) {
- case kMICRONS:
- par[14] = 0.1 * fThick * fDens * (fAmasr / fAmat);
- break;
- case kMGCM2:
- par[14] = fThick;
- break;
- case kCM:
- par[14] = fThick * fDens * 1000. * (fAmasr / fAmat);
- break;
- //gases: calculate density*thickness of cell
- case kMBAR:
- case kTORR:
- par[14] = fGasThick * GetDensity();
- break;
- }
- par[15] = A;
- par[16] = GetEminVedaloss(Z) * A;
- par[17] = GetEmaxVedaloss(Z) * A;
- if( fAmasr == fAmat ) par[18] = 0.;
- else par[18] = TMath::Log(fAmasr / fAmat);
-}
-
-//______________________________________________________________________________________________//
-void KVMaterial::SetELossParams(Int_t Z, Int_t A)
-{
- //Initialise energy loss coefficients for this material and a given incident nucleus (Z,A)
- static Double_t par[19];
- GetELossParams(Z, A, par);
- if (!ELoss) {
- //search in list of functions for one corresponding to this material
- //the name of the required function is ELoss
- ELoss = (TF1 *) gROOT->GetListOfFunctions()->FindObject("ELoss");
- if (!ELoss)
- ELoss = new TF1("ELoss", ELossSaclay, 0.1, 5000., 19);
- }
- ELoss->SetParameters(par);
-}
-
-//______________________________________________________________________________________________//
-void KVMaterial::SetEResParams(Int_t Z, Int_t A)
-{
- //Initialise coefficients for residual energy function for this material and a given incident nucleus (Z,A)
-
- static Double_t par[19];
- GetELossParams(Z, A, par);
- if (!ERes) {
- //search in list of functions for one corresponding to this material
- //the name of the required function is ERes
- ERes = (TF1 *) gROOT->GetListOfFunctions()->FindObject("ERes");
- if (!ERes)
- ERes = new TF1("ERes", EResSaclay, 0.1, 5000., 19);
- }
- ERes->SetParameters(par);
-}
-
//______________________________________________________________________________________//
Double_t KVMaterial::GetELostByParticle(KVNucleus * kvn, TVector3 * norm)
{
- //*** Saclay VEDALOSS library: translation of functions PARAM, GET_PARAM and ELOSS.
- //
- //Method to simulate passage of a charged particle (ion) through a given thickness of
- //material. The energy loss of the particle in the material is returned by the method.
- //The thickness is given in cm (for CsI crystals, NE102), in mbar (for C3F8), in micrometres
- //(for Silicon, Mylar), or mg/cm2 (all others).
- //
- //If the optional argument 'norm' is given, it is supposed to be a vector
- //normal to the material, oriented from the origin towards the material.
- //In this case the effective thickness of the material 'seen' by the particle
- //depending on its direction of motion is used for the calculation.
-
-
+ // Method to simulate passage of an ion through a given thickness of material.
+ // The energy loss of the particle (in MeV) in the material is returned by the method.
+ //
+ // If the optional argument 'norm' is given, it is supposed to be a vector
+ // normal to the material, oriented from the origin towards the material.
+ // In this case the effective thickness of the material 'seen' by the particle
+ // depending on its direction of motion is used for the calculation.
+
+ Double_t thickness;
if (norm) {
- Double_t thick = GetThickness();
TVector3 p = kvn->GetMomentum();
- Double_t e = GetEffectiveThickness((*norm), p);
-#ifdef DBG_TRGT
- cout << "Thickness used in energy loss calculation=" << e << endl;
-#endif
- SetThickness(e);
- Double_t E_loss = GetDeltaE(kvn->GetZ(), kvn->GetA(), kvn->GetKE());
- SetThickness(thick);
- return E_loss;
+ thickness = GetEffectiveThickness((*norm), p);
}
- Double_t E_loss = GetDeltaE(kvn->GetZ(), kvn->GetA(), kvn->GetKE());
+ else
+ thickness = GetThickness();
+ Double_t E_loss =
+ fIonRangeTable->GetLinearDeltaEOfIon(GetType(),kvn->GetZ(), kvn->GetA(), kvn->GetKE(),
+ thickness, fAmasr,fTemp,fPressure);
return E_loss;
}
@@ -804,16 +434,16 @@
//In this case the effective thickness of the material 'seen' by the particle
//depending on its direction of motion is used for the calculation.
+ Double_t thickness;
if (norm) {
- Double_t thick = GetThickness();
TVector3 p = kvn->GetMomentum();
- Double_t e = GetEffectiveThickness((*norm), p);
- SetThickness(e);
- Double_t E_inc = GetIncidentEnergyFromERes(kvn->GetZ(), kvn->GetA(), kvn->GetKE());
- SetThickness(thick);
- return E_inc;
+ thickness = GetEffectiveThickness((*norm), p);
}
- Double_t E_inc = GetIncidentEnergyFromERes(kvn->GetZ(), kvn->GetA(), kvn->GetKE());
+ else
+ thickness = GetThickness();
+ Double_t E_inc = fIonRangeTable->
+ GetLinearEIncFromEResOfIon(GetType(),kvn->GetZ(), kvn->GetA(), kvn->GetKE(),
+ thickness, fAmasr, fTemp, fPressure);
return E_inc;
}
@@ -823,14 +453,10 @@
{
// Calculate energy loss in absorber for incident nucleus (Z,A)
// with kinetic energy Einc (MeV)
- // For a KVDetector, the energy loss of the ACTIVE layer only is returned.
-
- if(Z<1 || Z>ZMAX_VEDALOSS) return 0.;
-
- SetELossParams(Z, A);
-
- Double_t E_loss = GetELossFunction()->Eval(Einc);
-
+
+ if(Z<1) return 0.;
+ Double_t E_loss =
+ fIonRangeTable->GetLinearDeltaEOfIon(GetType(), Z, A, Einc, GetThickness(), fAmasr, fTemp, fPressure);
return E_loss;
}
@@ -840,27 +466,17 @@
{
// Calculate energy loss in absorber for nucleus (Z,A)
// having a residual kinetic energy Eres (MeV) after the absorber
- // For a KVDetector, the energy loss of the ACTIVE layer only is returned.
-
- if(Z<1 || Z>ZMAX_VEDALOSS) return 0.;
-
- Double_t EINC = -1.;
- SetEResParams(Z, A);
- EINC = GetEResFunction()->GetX(Eres, GetEminVedaloss(Z) * A,
- GetEmaxVedaloss(Z) * A);
-
- // Calculate energy loss in absorber for incident nucleus (Z,A)
- // with kinetic energy Einc (MeV)
- SetELossParams(Z, A);
- Double_t E_loss = GetELossFunction()->Eval(EINC);
-
+
+ if(Z<1) return 0.;
+ Double_t E_loss = fIonRangeTable->
+ GetLinearDeltaEFromEResOfIon(
+ GetType(),Z,A,Eres,GetThickness(),fAmasr,fTemp,fPressure);
return E_loss;
}
//______________________________________________________________________________________//
-Double_t KVMaterial::GetEResFromDeltaE(Int_t Z, Int_t A, Double_t dE,
- enum SolType type)
+Double_t KVMaterial::GetEResFromDeltaE(Int_t Z, Int_t A, Double_t dE, enum KVIonRangeTable::SolType type)
{
// Calculate residual kinetic energy Eres (MeV) after the absorber from
// energy loss in absorber for nucleus (Z,A). If dE is not given, the energy
@@ -880,54 +496,29 @@
//maximum dE.
Double_t EINC = GetIncidentEnergy(Z,A,dE,type);
- SetEResParams(Z, A);
- return GetEResFunction()->Eval(EINC);
+ return EINC-dE;
}
//__________________________________________________________________________________________
-Double_t KVMaterial::GetIncidentEnergy(Int_t Z, Int_t A, Double_t delta_e,
- enum SolType type)
+Double_t KVMaterial::GetIncidentEnergy(Int_t Z, Int_t A, Double_t delta_e, enum KVIonRangeTable::SolType type)
{
//Calculate incident energy of nucleus (Z,A) corresponding to the energy loss
//in this absorber. If delta_e is not given, the energy loss in this absorber is used.
//
//By default the solution corresponding to the highest incident energy is returned
//This is the solution found for Einc greater than the maximum of the dE(Einc) curve.
- //If you want the low energy solution (i.e. below Bragg peak for gas detector)
- //set SolType = KVMaterial::kEmin.
+ //If you want the low energy solution set SolType = KVIonRangeTable::kEmin.
//
//If the given energy loss in the absorber is greater than the maximum theoretical dE
- //(dE > GetBraggDE) then we return the incident energy corresponding to the maximum,
- //GetBraggE.
+ //(dE > GetMaxDeltaE(Z,A)) then we return the incident energy corresponding to the maximum,
+ //GetEIncOfMaxDeltaE(Z,A)
- if(Z<1 || Z>ZMAX_VEDALOSS) return 0.;
+ if(Z<1) return 0.;
- Double_t EINC = -1.;
Double_t DE = (delta_e > 0 ? delta_e : GetEnergyLoss());
- SetELossParams(Z, A);
- //energy loss in the absorber is greater than the maximum theoretical dE ?
- Double_t braggDE = GetELossFunction()->GetMaximum(GetEminVedaloss(Z) * A,
- GetEmaxVedaloss(Z) * A);
- if (DE > braggDE)
- return GetELossFunction()->GetMaximumX(GetEminVedaloss(Z) * A,
- GetEmaxVedaloss(Z) * A);
-
- // Emax = GetBraggE(Z,A)
- Double_t Emax = GetELossFunction()->GetMaximumX(GetEminVedaloss(Z) * A,
- GetEmaxVedaloss(Z) * A);
- switch (type) {
- case kEmin:
- EINC = GetELossFunction()->GetX(DE, GetEminVedaloss(Z) * A, Emax);
- break;
-
- case kEmax:
- EINC = GetELossFunction()->GetX(DE, Emax, GetEmaxVedaloss(Z) * A);
- break;
- }
-
- return EINC;
+ return fIonRangeTable->GetLinearEIncFromDeltaEOfIon(GetType(),Z,A,DE,GetThickness(),type,fAmasr,fTemp,fPressure);
}
//______________________________________________________________________________________//
@@ -937,11 +528,10 @@
// Calculate residual energy after absorber for incident nucleus (Z,A)
// with kinetic energy Einc (MeV)
- if(Z<1 || Z>ZMAX_VEDALOSS) return 0.;
+ if(Z<1) return 0.;
- SetEResParams(Z, A);
-
- Double_t E_res = GetEResFunction()->Eval(Einc);
+ Double_t E_res =
+ fIonRangeTable->GetEResOfIon(GetType(), Z, A, Einc, fThick, fAmasr, fTemp, fPressure);
return E_res;
}
@@ -985,20 +575,6 @@
fELoss = 0.0;
}
-//__________________________________________________________________________________________
-
-const Char_t *KVMaterial::GetThicknessUnits() const
-{
- //Returns a string with the name of the units used to measure "thickness"
- //for this material type
- if (GetActiveLayer())
- return GetActiveLayer()->GetThicknessUnits();
- if (fUnits < 10)
- return kUnits[GetUnits()];
- else
- return "(undefined)";
-}
-
#if ROOT_VERSION_CODE >= ROOT_VERSION(3,4,0)
void KVMaterial::Copy(TObject & obj) const
#else
@@ -1009,9 +585,9 @@
KVBase::Copy(obj);
((KVMaterial &) obj).SetMaterial(GetType());
((KVMaterial &) obj).SetMass(GetMass());
- ((KVMaterial &) obj).SetThickness(GetThickness());
((KVMaterial &) obj).SetPressure(GetPressure());
((KVMaterial &) obj).SetTemperature(GetTemperature());
+ ((KVMaterial &) obj).SetThickness(GetThickness());
}
//__________________________________________________________________________________________
@@ -1020,39 +596,36 @@
Double_t Eres)
{
//Get incident energy from residual energy
- if(Z<1 || Z>ZMAX_VEDALOSS) return 0.;
+ if(Z<1) return 0.;
- Double_t EINC = -1.;
- SetEResParams(Z, A);
- EINC =
- GetEResFunction()->GetX(Eres, GetEminVedaloss(Z) * A,
- GetEmaxVedaloss(Z) * A);
- return EINC;
+ return fIonRangeTable->
+ GetLinearEIncFromEResOfIon(GetType(),Z,A,Eres,GetThickness(),fAmasr,fTemp,fPressure);
}
//__________________________________________________________________________________________
-Double_t KVMaterial::GetBraggE(Int_t Z, Int_t A)
+Double_t KVMaterial::GetEIncOfMaxDeltaE(Int_t Z, Int_t A)
{
//Incident energy for which the DE(E) curve has a maximum
- if(Z<1 || Z>ZMAX_VEDALOSS) return 0.;
+ if(Z<1) return 0.;
- SetELossParams(Z, A);
- return GetELossFunction()->GetMaximumX(GetEminVedaloss(Z) * A,
- GetEmaxVedaloss(Z) * A);
+ return fIonRangeTable->
+ GetLinearEIncOfMaxDeltaEOfIon(GetType(),Z,A,GetThickness(),fAmasr,fTemp,fPressure);
}
//__________________________________________________________________________________________
-Double_t KVMaterial::GetBraggDE(Int_t Z, Int_t A)
+Double_t KVMaterial::GetMaxDeltaE(Int_t Z, Int_t A)
{
- //The maximum energy loss of this particle
- //(corresponding to incident energy GetBraggE(Z,A))
- if(Z<1 || Z>ZMAX_VEDALOSS) return 0.;
-
- SetELossParams(Z, A);
- return GetELossFunction()->GetMaximum(GetEminVedaloss(Z) * A,
- GetEmaxVedaloss(Z) * A);
+ // The maximum energy loss of this particle (corresponding to incident energy GetEIncOfMaxDeltaE(Z,A))
+ // For detectors, this is the maximum energy loss in the active layer.
+
+ if(GetActiveLayer()) return GetActiveLayer()->GetMaxDeltaE(Z,A);
+
+ if(Z<1) return 0.;
+
+ return fIonRangeTable->
+ GetLinearMaxDeltaEOfIon(GetType(),Z,A,GetThickness(),fAmasr,fTemp,fPressure);
}
//__________________________________________________________________________________________
@@ -1108,30 +681,12 @@
return gmed;
}
-
-//______________________________________________________________________________
-
-void KVMaterial::Streamer(TBuffer &R__b)
+Double_t KVMaterial::GetEmaxValid(Int_t Z, Int_t A)
{
- // Stream an object of class KVMaterial.
+ // Return maximum incident energy for which range tables are valid
+ // for this material and ion (Z,A).
+ // For detectors, return max energy for active layer.
+ if(GetActiveLayer()) return GetActiveLayer()->GetEmaxValid(Z,A);
+ return fIonRangeTable->GetEmaxValid(GetType(), Z, A);
+}
- UInt_t R__s, R__c;
- if (R__b.IsReading()) {
- Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
- R__b.ReadClassBuffer(KVMaterial::Class(),this,R__v,R__s,R__c);
- if (R__v < 5){
- // correct densities of solids written with class version < 5
- // this is to fix launchpad bug#446163
- // it affects materials with fUnits = kMICRONS or kCM
- // if the density we read is less than that we find
- // for a new instanciation of the same type of material,
- // we use the current value.
- if(fUnits==kMICRONS || fUnits==kCM){
- KVMaterial tmp(GetType());
- if(GetDensity() < tmp.GetDensity()) SetDensity(tmp.GetDensity());
- }
- }
- } else {
- R__b.WriteClassBuffer(KVMaterial::Class(),this);
- }
-}
=== modified file 'KVMultiDet/detectors/KVMaterial.h'
--- KVMultiDet/detectors/KVMaterial.h 2009-11-26 16:48:10 +0000
+++ KVMultiDet/detectors/KVMaterial.h 2011-03-14 15:58:00 +0000
@@ -1,5 +1,4 @@
/***************************************************************************
-$Id: KVMaterial.h,v 1.31 2008/12/17 13:01:26 franklan Exp $
kvmaterial.h - description
-------------------
begin : Thu May 16 2002
@@ -22,115 +21,52 @@
#define KVMATERIAL_H
#include "RVersion.h"
-#include "TRef.h"
#include "KVBase.h"
#include "KVList.h"
#include "TF1.h"
#include "TVector3.h"
#include "Riostream.h"
-
-// warning messages
-#define KVMATERIAL_INIT_FILE_NOT_FOUND "Initialisation file %s not found. Cannot create detectors"
-#define KVMATERIAL_INIT_FILE_HEADER_PROBLEM "Problem reading initialisation file %s. Check format of header"
-#define KVMATERIAL_INIT_FILE_COEFF_PROBLEM "Problem reading initialisation file %s. Check format of coefficients for %s"
-#define KVMATERIAL_INIT_FILE_DATA_PROBLEM "Problem reading initialisation file %s. Check format of data"
-#define KVMATERIAL_TYPE_NOT_FOUND "Unknown material type %s. Check file %s for known types"
-
-#define RTT 623.61040
-#define ZERO_KELVIN 273.15
-#define PERC 0.02
-
-// maximum atomic number included in range tables
-#define ZMAX_VEDALOSS 100
+#include "KVIonRangeTable.h"
class KVNucleus;
class KVTelescope;
class TGeoMedium;
-Double_t ELossSaclay(Double_t * x, Double_t * par);
-Double_t EResSaclay(Double_t * x, Double_t * par);
-
class KVMaterial:public KVBase {
protected:
- TF1 * ELoss; //!function calculating energy loss in material from incident energy
- TF1 *ERes; //!function calculating residual energy from incident energy
+ static KVIonRangeTable* fIonRangeTable; // pointer to class used to calculate charged particle ranges & energy losses
+
TVector3 fNormToMat;//!dummy vector for calculating normal to absorber
public:
- TF1 * GetELossFunction() {
- return ELoss;
- };
- TF1 *GetEResFunction() {
- return ERes;
- };
- virtual void SetELossParams(Int_t Z, Int_t A);
- void GetELossParams(Int_t Z, Int_t A, Double_t * par);
- enum SolType {
- kEmax,
- kEmin
- };
- virtual void SetEResParams(Int_t Z, Int_t A);
-
- virtual Double_t GetIncidentEnergyFromERes(Int_t Z, Int_t A,
- Double_t Eres);
-
- enum { //constants defined for units
- kMBAR,
- kMICRONS,
- kMGCM2,
- kCM,
- kTORR
- };
private:
- Double_t fAmat; //atomic mass of material
- Double_t fZmat; //atomic number of material
- Double_t fDens; //density of material
- UInt_t fUnits; //used to determine how to calculate effective thickness
- Double_t fGasThick; //only used for gases
- Double_t fMoleWt; //only used for gases
- Double_t fTemp; //only used for gases
- Double_t fCoeff[ZMAX_VEDALOSS][15]; //[ZMAX_VEDALOSS][15] parameters for range tables
- Float_t fAmasr; //isotopic mass for material; if not set will use default value
- static TString kDataFilePath; //! full path to file containing energy loss parameters
- Float_t fEmax[ZMAX_VEDALOSS]; //[ZMAX_VEDALOSS] Z-dependent maximum energy/nucleon for calculation to be valid
- Float_t fEmin[ZMAX_VEDALOSS]; //[ZMAX_VEDALOSS] Z-dependent minimum energy/nucleon for calculation to be valid
- Float_t fThick; //thickness of absorber
+ Int_t fAmasr; // isotopic mass of element
+ Double_t fThick; // area density of absorber in g/cm**2
+ Double_t fPressure; // gas pressure in torr
+ Double_t fTemp; // gas temperature in degrees celsius
Double_t fELoss; //total of energy lost by all particles traversing absorber
- static Char_t kUnits[5][10]; //!strings for names of units for thickness
public:
KVMaterial();
- KVMaterial(const Char_t * type, const Float_t thick = 0.0);
+ KVMaterial(const Char_t * type, const Double_t thick = 0.0);
+ KVMaterial(const Char_t * gas, const Double_t thick, const Double_t pressure, const Double_t temperature = 19.0);
+ KVMaterial(Double_t area_density, const Char_t * type);
KVMaterial(const KVMaterial &);
void init();
virtual ~ KVMaterial();
- static KVList *GetListOfMaterials();
- void SetMass(Float_t a);
+ void SetMass(Double_t a);
virtual void SetMaterial(const Char_t * type);
Double_t GetMass() const;
Double_t GetZ() const;
Double_t GetDensity() const;
- void SetDensity(Double_t d)
- {
- // Set density of material. Units are g/cm**3
- fDens = d;
- }
- Float_t GetEmaxVedaloss(UInt_t Z) const {
- if (GetActiveLayer())
- return GetActiveLayer()->GetEmaxVedaloss(Z);
- return ((Z <= ZMAX_VEDALOSS && Z > 0) ? fEmax[Z - 1] : 0.0);
- };
- Float_t GetEminVedaloss(UInt_t Z) const {
- if (GetActiveLayer())
- return GetActiveLayer()->GetEminVedaloss(Z);
- return ((Z <= ZMAX_VEDALOSS && Z > 0) ? fEmin[Z - 1] : 0.0);
- };
- void SetThickness(Float_t thick);
- Float_t GetThickness() const;
- Double_t GetThicknessInCM() const;
+ void SetAreaDensity(Double_t dens /* g/cm**2 */);
+ Double_t GetAreaDensity() const;
+ virtual void SetThickness(Double_t thick /* cm */);
+ virtual Double_t GetThickness() const;
Double_t GetEffectiveThickness(TVector3 & norm, TVector3 & direction);
+ Double_t GetEffectiveAreaDensity(TVector3 & norm, TVector3 & direction);
virtual void DetectParticle(KVNucleus *, TVector3 * norm = 0);
virtual Double_t GetELostByParticle(KVNucleus *, TVector3 * norm = 0);
@@ -151,33 +87,33 @@
#endif
virtual void Clear(Option_t * opt = "");
- const Char_t *GetThicknessUnits() const;
- virtual UInt_t GetUnits() const;
- void SetUnits(UInt_t);
-
- Double_t GetIncidentEnergy(Int_t Z, Int_t A, Double_t delta_e =
- -1.0, enum SolType type = kEmax);
- Double_t GetDeltaE(Int_t Z, Int_t A, Double_t Einc);
- Double_t GetDeltaEFromERes(Int_t Z, Int_t A, Double_t Eres);
- Double_t GetERes(Int_t Z, Int_t A, Double_t Einc);
- Double_t GetEResFromDeltaE(Int_t Z, Int_t A, Double_t dE =
- -1.0, enum SolType type = kEmax);
- Double_t GetBraggE(Int_t Z, Int_t A);
- Double_t GetBraggDE(Int_t Z, Int_t A);
+ virtual Double_t GetEmaxValid(Int_t Z, Int_t A);
+ virtual Double_t GetIncidentEnergy(Int_t Z, Int_t A, Double_t delta_e =
+ -1.0, enum KVIonRangeTable::SolType type = KVIonRangeTable::kEmax);
+ virtual Double_t GetIncidentEnergyFromERes(Int_t Z, Int_t A, Double_t Eres);
+ virtual Double_t GetDeltaE(Int_t Z, Int_t A, Double_t Einc);
+ virtual Double_t GetDeltaEFromERes(Int_t Z, Int_t A, Double_t Eres);
+ virtual Double_t GetERes(Int_t Z, Int_t A, Double_t Einc);
+ virtual Double_t GetEResFromDeltaE(Int_t Z, Int_t A, Double_t dE =
+ -1.0, enum KVIonRangeTable::SolType type = KVIonRangeTable::kEmax);
+ virtual Double_t GetEIncOfMaxDeltaE(Int_t Z, Int_t A);
+ virtual Double_t GetMaxDeltaE(Int_t Z, Int_t A);
virtual KVMaterial *GetActiveLayer() const {
return 0;
}; //to facilitate polymorphism with KVDetector class
- virtual void SetPressure(Float_t);
- virtual void SetTemperature(Float_t t);
+ virtual void SetPressure(Double_t);
+ virtual void SetTemperature(Double_t);
- Float_t GetPressure() const;
- Float_t GetTemperature() const;
+ virtual Double_t GetPressure() const;
+ virtual Double_t GetTemperature() const;
Bool_t IsIsotopic() const;
Bool_t IsNat() const;
+ Bool_t IsGas() const;
+
virtual const TVector3& GetNormal() {
// Return vector normal to surface of absorber. For a KVMaterial, this is (0,0,1) as
// a basic absorber has no orientation. Rederived in child classes KVTarget and KVDetector.
@@ -185,8 +121,10 @@
};
virtual TGeoMedium* GetGeoMedium(const Char_t* /*med_name*/="");
+
+ virtual KVIonRangeTable* GetRangeTable() const { return fIonRangeTable; };
- ClassDef(KVMaterial, 5) //Interface to VEDALOSS description of slowing-down of ions in absorbers
+ ClassDef(KVMaterial, 6) // Class describing physical materials used to construct detectors & targets
};
#endif
=== modified file 'KVMultiDet/detectors/KVTarget.cpp'
--- KVMultiDet/detectors/KVTarget.cpp 2009-11-29 18:21:06 +0000
+++ KVMultiDet/detectors/KVTarget.cpp 2011-03-14 15:58:00 +0000
@@ -28,6 +28,9 @@
KVTarget targ("Ta", 8.3); // first layer 8.3mg/cm2 Ta
targ.AddLayer("C", 20.0); // second layer 20mg/cm2 C
</pre>
+
+Note that all "thicknesses" for targets and target layers are in mg/cm**2.
+
<h4>ORIENTATION OF TARGET</h4>
The target can be oriented in an arbitrary way, by defining the vector normal to its
@@ -139,10 +142,11 @@
//________________________________________________________________________
-KVTarget::KVTarget(const Char_t * material,
- const Double_t thick):KVMaterial(), fNormal(0, 0, 1)
+KVTarget::KVTarget(const Char_t * material, const Double_t thick):KVMaterial(), fNormal(0, 0, 1)
{
- //Just give the type & thickness of material for target
+ // Just give the type & "thickness" of material for target
+ // The "thickness" is the area density of the target in mg/cm**2.
+
init();
AddLayer(material, thick);
}
@@ -176,19 +180,19 @@
TIter next(GetLayers());
KVMaterial *mat = 0;
while ((mat = (KVMaterial *) next())) {
- ((KVTarget &) obj).AddLayer(mat->GetType(),mat->GetThickness());
+ ((KVTarget &) obj).AddLayer(mat->GetType(), mat->GetAreaDensity()/KVUnits::mg);
}
}
//________________________________________________________________________
void KVTarget::AddLayer(const Char_t * material, Double_t thick)
{
- //Add a layer to a target
- //Sets/updates name of target with name of material.
- //In case of multi-layer target the name is
+ // Add a layer to a target, with 'thickness' in mg/cm**2 (area density).
+ // Sets/updates name of target with name of material.
+ // In case of multi-layer target the name is
// material1/material2/material3/...
- fTargets->Add(new KVMaterial(material, thick));
+ fTargets->Add(new KVMaterial(thick*KVUnits::mg, material));
fNLayers++;
if (fNLayers == 1) {
SetName(material);
@@ -204,21 +208,25 @@
Double_t KVTarget::GetTotalThickness()
{
- //return sum of thicknesses of all layers in target
+ // return sum of 'thicknesses' (area densities in mg/cm**2)
+ // of all layers in target
+
Float_t thick = 0.;
TIter next(fTargets);
KVMaterial *mat = 0;
while ((mat = (KVMaterial *) next())) {
- thick += mat->GetThickness();
+ thick += mat->GetAreaDensity();
}
- return thick;
+ return thick/KVUnits::mg;
}
//________________________________________________________________________
Double_t KVTarget::GetTotalThickness(Int_t lay1, Int_t lay2)
{
- //return sum of thicknesses of layers lay1 to lay2 in target
+ //return sum of 'thicknesses' (area densities in mg/cm**2)
+ // of layers lay1 to lay2 in target
+
Double_t thick = 0.;
for (register int i = lay1; i <= lay2; i++) {
thick += GetThickness(i);
@@ -250,7 +258,7 @@
Double_t KVTarget::GetEffectiveThickness(KVParticle * part, Int_t ilayer)
{
- //Return effective thickness of layer ilayer (ilayer=1, 2, ...)
+ //Return effective 'thickness' (in mg/cm**2) of layer ilayer (ilayer=1, 2, ...)
//By default ilayer=1 (i.e. for single layer target)
//The effective thickness depends on the angle of the target (rotation about
//x-axis => theta wrt z- (beam)-axis).
@@ -273,7 +281,7 @@
Double_t KVTarget::GetEffectiveThickness(TVector3 & direction,
Int_t ilayer)
{
- //Return effective thickness of layer ilayer (ilayer=1, 2, ...)
+ //Return effective 'thickness' (in mg/cm**2) of layer ilayer (ilayer=1, 2, ...)
//By default ilayer=1 (i.e. for single layer target)
//The effective thickness depends on the orientation of the target (given by
//the direction of the normal to its surface) and on the direction (e.g. direction of a particle)
@@ -284,8 +292,7 @@
ilayer, NumberOfLayers());
return 0.0;
}
- return GetLayerByIndex(ilayer)->GetEffectiveThickness(fNormal,
- direction);
+ return GetLayerByIndex(ilayer)->GetEffectiveAreaDensity(fNormal, direction)/KVUnits::mg;
}
//________________________________________________________________________
@@ -294,7 +301,7 @@
{
//Returns absorber corresponding to 'depth' inside target, starting from the 'entrance'
//layer and following the direction of 'depth'. Note: 'depth' is measured in the same
- //'thickness' units as the thickness of the different layers of the target (mg/cm2, um, etc.)
+ //'thickness' units as the thickness of the different layers of the target (mg/cm2)
//WARNING : returns 0 if no layer is found (depth is outside of target)
return GetLayerByIndex(GetLayerIndex(depth));
@@ -306,7 +313,7 @@
{
//Returns absorber index corresponding to 'depth' inside target, starting from the 'entrance'
//layer and following the direction of 'depth'. Note: 'depth' is measured in the same
- //'thickness' units as the thickness of the different layers of the target (mg/cm2, um, etc.)
+ //'thickness' units as the thickness of the different layers of the target (mg/cm2)
//WARNING : user should check returned index is >0
//If not, this means that the given depth does not correspond to a layer inside the target
@@ -357,9 +364,9 @@
Double_t KVTarget::GetThickness(Int_t ilayer) const
{
- //Thickness of layer 'ilayer' in target
+ //'Thickness' in mg/cm**2 of layer 'ilayer' in target
KVMaterial *lay = GetLayerByIndex(ilayer);
- return (lay ? lay->GetThickness() : 0.0);
+ return (lay ? lay->GetAreaDensity()/KVUnits::mg : 0.0);
}
//________________________________________________________________________
@@ -368,7 +375,7 @@
{
//Returns absorber index corresponding to 'depth' inside target, starting from the 'entrance'
//layer and following the normal direction. Note: 'depth' is measured in the same
- //'thickness' units as the thickness of the different layers of the target (mg/cm2, um, etc.)
+ //'thickness' units as the thickness of the different layers of the target (mg/cm2)
//WARNING : user should check returned index is >0
//If not, this means that the given depth does not correspond to a layer inside the target
@@ -387,7 +394,7 @@
Double_t KVTarget::GetTotalEffectiveThickness(KVParticle * part)
{
- //return sum of effective thicknesses of all layers in target
+ //return sum of effective 'thicknesses' (mg/cm**2) of all layers in target
//taking into account the angle of the target to the beam
//and the direction of motion of the incident particle.
//If no particle is given, effective thicknesses are calculated as for
@@ -402,7 +409,7 @@
Double_t KVTarget::GetTotalEffectiveThickness(TVector3 & dir, Int_t ilay1,
Int_t ilay2)
{
- //return sum of effective thicknesses of layers ilay1 to ilay2 in target
+ //return sum of effective 'thicknesses' (mg/cm**2) of layers ilay1 to ilay2 in target
//taking into account the angle of the target to the beam
//and the given direction.
//
@@ -484,18 +491,18 @@
iplay_index) -
GetInteractionPoint() * GetNormal();
e_thick_iplay =
- (IsIncoming()? iplay->GetThickness() -
+ (IsIncoming()? iplay->GetAreaDensity()/KVUnits::mg -
e_thick_iplay : e_thick_iplay);
if (backwards)
- e_thick_iplay = iplay->GetThickness() - e_thick_iplay;
+ e_thick_iplay = iplay->GetAreaDensity()/KVUnits::mg - e_thick_iplay;
#ifdef DBG_TRGT
cout << "Effective thickness of IP layer is " << e_thick_iplay <<
- " (real:" << iplay->GetThickness() << ")" << endl;
+ " (real:" << iplay->GetAreaDensity()/KVUnits::mg << ")" << endl;
#endif
//modify effective physical thickness of layer
- Double_t thick_iplay = iplay->GetThickness();
- iplay->SetThickness(e_thick_iplay);
+ Double_t thick_iplay = iplay->GetAreaDensity();// in g/cm**2
+ iplay->SetAreaDensity(e_thick_iplay*KVUnits::mg);
//first and last indices of layers to pass through
Int_t ilay1 =
@@ -516,7 +523,7 @@
}
//reset original thickness of IP layer
- iplay->SetThickness(thick_iplay);
+ iplay->SetAreaDensity(thick_iplay);
} else {
Error("DetectParticle", "Interaction point is outside of target");
@@ -586,18 +593,18 @@
iplay_index) -
GetInteractionPoint() * GetNormal();
e_thick_iplay =
- (IsIncoming()? iplay->GetThickness() -
+ (IsIncoming()? iplay->GetAreaDensity()/KVUnits::mg -
e_thick_iplay : e_thick_iplay);
if (backwards)
- e_thick_iplay = iplay->GetThickness() - e_thick_iplay;
+ e_thick_iplay = iplay->GetAreaDensity()/KVUnits::mg - e_thick_iplay;
#ifdef DBG_TRGT
cout << "Effective thickness of IP layer is " << e_thick_iplay <<
- " (real:" << iplay->GetThickness() << ")" << endl;
+ " (real:" << iplay->GetAreaDensity()/KVUnits::mg << ")" << endl;
#endif
//modify effective physical thickness of layer
- Double_t thick_iplay = iplay->GetThickness();
- iplay->SetThickness(e_thick_iplay);
+ Double_t thick_iplay = iplay->GetAreaDensity(); // g/cm**2
+ iplay->SetAreaDensity(e_thick_iplay*KVUnits::mg);
//first and last indices of layers to pass through
Int_t ilay1 =
@@ -626,7 +633,7 @@
}
//reset original thickness of IP layer
- iplay->SetThickness(thick_iplay);
+ iplay->SetAreaDensity(thick_iplay);
} else {
Error("DetectParticle", "Interaction point is outside of target");
@@ -776,18 +783,18 @@
void KVTarget::SetMaterial(const Char_t * type)
{
- //Set material of first layer
+ // Set material of first layer
if (GetLayerByIndex(1))
GetLayerByIndex(1)->SetMaterial(type);
}
//______________________________________________________________________________________________________
-void KVTarget::SetThickness(Float_t thick, Int_t ilayer)
+void KVTarget::SetLayerThickness(Float_t thick, Int_t ilayer)
{
- //Setthickness of a layer, by default this is the first layer
+ // Set 'thickness' in mg/cm**2 of a layer, by default this is the first layer
if (GetLayerByIndex(ilayer))
- GetLayerByIndex(ilayer)->SetThickness(thick);
+ GetLayerByIndex(ilayer)->SetAreaDensity(thick*KVUnits::mg);
}
//______________________________________________________________________________________________________
@@ -795,7 +802,7 @@
Double_t KVTarget::GetAtomsPerCM2() const
{
//Calculates total number of atoms per square centimetre of the target.
- //For a multilayer target, the areal densities for each layer are summed up.
+ //For a multilayer target, the area densities for each layer are summed up.
Double_t atom_cib = 0;
for (register int i = 1; i <= NumberOfLayers(); i++) {
//N_atoms = N_Avogadro * target_thickness (mg/cm**2) * 1.e-3 / atomic_mass_of_target
@@ -860,14 +867,14 @@
iplay_index) -
GetInteractionPoint() * GetNormal();
e_thick_iplay =
- (IsIncoming()? iplay->GetThickness() -
+ (IsIncoming()? iplay->GetAreaDensity()/KVUnits::mg -
e_thick_iplay : e_thick_iplay);
if (backwards)
- e_thick_iplay = iplay->GetThickness() - e_thick_iplay;
+ e_thick_iplay = iplay->GetAreaDensity()/KVUnits::mg - e_thick_iplay;
//modify effective physical thickness of layer
Double_t thick_iplay = iplay->GetThickness();
- iplay->SetThickness(e_thick_iplay);
+ iplay->SetAreaDensity(e_thick_iplay);
//first and last indices of layers to pass through
Int_t ilay1 =
@@ -901,7 +908,7 @@
}
//reset original thickness of IP layer
- iplay->SetThickness(thick_iplay);
+ iplay->SetAreaDensity(thick_iplay);
} else {
Error("GetParticleEIncFromERes", "Interaction point is outside of target");
@@ -940,20 +947,3 @@
return GetParticleEIncFromERes(&dummy);
}
-//___________________________________________________________________________________
-
-UInt_t KVTarget::GetUnits() const
-{
- // Returns units for 'thickness' of target.
- // We actually only look at the 1st layer of the target, and assume all layers have
- // the same units.
- // If you want, you can examine individual layers with GetLayer(i)->GetUnits().
- // Returned value can be one of:
- // KVMaterial::kMBAR, KVMaterial::kMGCM2, KVMaterial::kTORR,
- // KVMaterial::kMICRON, KVMaterial::kCM.
-
- KVMaterial * lay = const_cast<KVTarget*>(this)->GetLayerByIndex(1);
- if (lay) return lay->GetUnits();
- return 10; // no layers - undefined
-}
-
=== modified file 'KVMultiDet/detectors/KVTarget.h'
--- KVMultiDet/detectors/KVTarget.h 2009-11-26 16:48:10 +0000
+++ KVMultiDet/detectors/KVTarget.h 2011-03-14 15:58:00 +0000
@@ -49,7 +49,7 @@
Double_t GetAngleToBeam();
virtual void SetMaterial(const Char_t * type);
- void SetThickness(Float_t thick, Int_t ilayer = 1);
+ void SetLayerThickness(Float_t thick, Int_t ilayer = 1);
void AddLayer(const Char_t * material, Double_t thick);
Int_t NumberOfLayers() const {
@@ -73,9 +73,8 @@
Double_t GetTotalThickness();
Double_t GetTotalThickness(Int_t lay1, Int_t lay2);
- Float_t GetThickness() const {
- return (Float_t) const_cast <
- KVTarget * >(this)->GetTotalThickness();
+ Double_t GetThickness() const {
+ return const_cast <KVTarget * >(this)->GetTotalThickness();
};
Double_t GetTotalEffectiveThickness(KVParticle * part = 0);
Double_t GetTotalEffectiveThickness(TVector3 &, Int_t lay1 =
@@ -149,7 +148,7 @@
#else
virtual void Copy(TObject & obj);
#endif
- virtual UInt_t GetUnits() const;
+ // virtual UInt_t GetUnits() const;
Double_t GetAtomsPerCM2() const;
=== modified file 'KVMultiDet/events/KVElasticScatterEvent.cpp'
--- KVMultiDet/events/KVElasticScatterEvent.cpp 2011-02-17 07:50:39 +0000
+++ KVMultiDet/events/KVElasticScatterEvent.cpp 2011-03-14 15:58:00 +0000
@@ -401,10 +401,9 @@
);
printf("# Propagation dans une cible de:\n");
for (Int_t nn=0;nn<GetTarget().GetLayers()->GetEntries();nn+=1){
- printf("# type:%s epaisseur:%1.2lf (%s)\n",
+ printf("# type:%s epaisseur:%1.2lf (mg/cm**2)\n",
GetTarget().GetLayerByIndex(nn+1)->GetType(),
- GetTarget().GetLayerByIndex(nn+1)->GetThickness(),
- GetTarget().GetLayerByIndex(nn+1)->GetThicknessUnits()
+ GetTarget().GetLayerByIndex(nn+1)->GetAreaDensity()/(KVUnits::mg/pow(KVUnits::cm, 2.))
);
}
printf("#####################\n");
=== modified file 'KVMultiDet/gui/KVDBSystemDialog.cpp'
--- KVMultiDet/gui/KVDBSystemDialog.cpp 2010-01-07 14:25:47 +0000
+++ KVMultiDet/gui/KVDBSystemDialog.cpp 2011-03-14 15:58:00 +0000
@@ -435,12 +435,12 @@
return;
}
fLayer = (KVMaterial*)fTarget->GetLayers()->At(ind);
- //update thickness
- fNumberEntry1526->SetNumber(fLayer->GetThickness());
+ //update thickness - actually area density in mg/cm2
+ fNumberEntry1526->SetNumber(fLayer->GetAreaDensity()/(KVUnits::mg/pow(KVUnits::cm,2)));
//update atomic mass
fNumberEntry1537->SetNumber(fLayer->GetMass());
//update thickness units
- fLabel1530->SetText(fLayer->GetThicknessUnits());
+ fLabel1530->SetText("mg/cm2");
//enable button to remove layer
fTextButton1554->SetEnabled(kTRUE);
}
@@ -505,10 +505,11 @@
void KVDBSystemDialog::TargetLayerThicknessChanged(Long_t)
{
- //Called when target layer thickness is changed
+ // Called when target layer "thickness" is changed
+ // Note that this is in fact the area density in mg/cm**2
Double_t t = fNumberEntry1526->GetNumber();
- fLayer->SetThickness(t);
+ fLayer->SetAreaDensity(t*KVUnits::mg/pow(KVUnits::cm,2.));
SetNeedSave(1);
}
@@ -538,7 +539,9 @@
//The new layer will be added after any existing layers in the target.
//get selected material
- KVMaterial* mat = (KVMaterial*)fMaterialsList->At( fComboBox1542->GetSelected() );
+ TNamed* mat = (TNamed*)fMaterialsList->At( fComboBox1542->GetSelected() );
+ KVMaterial bidon;
+ KVIonRangeTable* RT = bidon.GetRangeTable();
//add to target of current system
//if no target is defined, we create a new one
if( !fTarget ){
@@ -546,16 +549,16 @@
if(fSystem){
fSystem->SetTarget( fTarget );
cout << "Created target for system : " << fSystem->GetName() << endl;
- fSystem->SetZtarget((UInt_t)mat->GetZ()); fSystem->SetAtarget((UInt_t)mat->GetMass());
+ fSystem->SetZtarget((UInt_t)RT->GetZ(mat->GetName())); fSystem->SetAtarget((UInt_t)RT->GetAtomicMass(mat->GetName()));
}
}
- //add layer with default thickness 0.1 mg/cm2
- fTarget->AddLayer( mat->GetType() , 0.1 );
+ //add layer with default area density 0.1 mg/cm2
+ fTarget->AddLayer( mat->GetTitle() , 0.1*KVUnits::mg/pow(KVUnits::cm,2.) );
//update list of layers in target
Int_t nlay = fComboBox1515->GetNumberOfEntries();
fComboBox1515->AddEntry( mat->GetName(), nlay );
fComboBox1515->Select( nlay+1 );
- cout << "Added layer " << mat->GetType() << " to target" << endl;
+ cout << "Added layer " << mat->GetTitle() << " to target" << endl;
SetNeedSave(1);
}
@@ -575,7 +578,7 @@
TIter next(fTarget->GetLayers()); KVMaterial* mat;
while ( (mat = (KVMaterial*)next()) ) {
if( mat != fLayer ){
- new_target->AddLayer( mat->GetType(), mat->GetThickness() );
+ new_target->AddLayer( mat->GetType(), mat->GetAreaDensity() );
}
}
new_target->SetAngleToBeam( fTarget->GetAngleToBeam() );
@@ -799,7 +802,8 @@
fComboBox1542 = new TGComboBox(fGroupFrame1541,-1,kHorizontalFrame | kSunkenFrame | kDoubleBorder | kOwnBackground);
fComboBox1542->Resize(80,22);
//fill list of all available materials
- fMaterialsList = KVMaterial::GetListOfMaterials();
+ KVMaterial bidon;
+ fMaterialsList = bidon.GetRangeTable()->GetListOfMaterials();
TIter it_mat(fMaterialsList); TObject* obj; Int_t ind = 0;
while( (obj = it_mat()) ){
fComboBox1542->AddEntry(obj->GetName(), ind++);
=== modified file 'KVMultiDet/gui/KVDBSystemDialog.h'
--- KVMultiDet/gui/KVDBSystemDialog.h 2009-01-22 15:03:32 +0000
+++ KVMultiDet/gui/KVDBSystemDialog.h 2011-03-14 15:58:00 +0000
@@ -57,7 +57,7 @@
KVMaterial *fLayer;//current target layer
KVNumberList fRuns;//runs selected by user
- KVList* fMaterialsList;//list of all available materials
+ TObjArray* fMaterialsList;//list of all available materials
KVDatedFileManager* fUndo;//allows to undo changes to Systems.dat
KVString fCurrentSystemsFile;//name (including timestamp) of currently used Systems.dat
=== modified file 'KVMultiDet/gui/KVDetectorBrowser.cpp'
--- KVMultiDet/gui/KVDetectorBrowser.cpp 2009-01-22 15:03:32 +0000
+++ KVMultiDet/gui/KVDetectorBrowser.cpp 2011-03-14 15:58:00 +0000
@@ -29,11 +29,11 @@
widg->SetWidget(fCombo);
AddToWidgetList(widg);
// fill list with all available materials
- KVList *mat_list = fDet->GetListOfMaterials();
+ TObjArray *mat_list = fDet->GetRangeTable()->GetListOfMaterials();
TIter next_mat(mat_list);
- KVMaterial *mat;
+ TNamed *mat;
Int_t i = 0, iselect = -1;
- while ((mat = (KVMaterial *) next_mat())) {
+ while ((mat = (TNamed *) next_mat())) {
if (!strcmp(mat->GetName(), fDet->GetMaterialName()))
iselect = i;
fCombo->AddEntry(mat->GetName(), i++);
=== modified file 'KVMultiDet/idmaps/KVIDZAGrid.cpp'
--- KVMultiDet/idmaps/KVIDZAGrid.cpp 2011-02-23 16:30:57 +0000
+++ KVMultiDet/idmaps/KVIDZAGrid.cpp 2011-03-14 15:58:00 +0000
@@ -1464,8 +1464,8 @@
Double_t E1, E2;
//find E1
- //go from 0.1 MeV to dE->GetBraggE(part.GetZ(),part.GetA()))
- Double_t E1min = 0.1, E1max = dEDet->GetBraggE(part.GetZ(),part.GetA());
+ //go from 0.1 MeV to dE->GetEIncOfMaxDeltaE(part.GetZ(),part.GetA()))
+ Double_t E1min = 0.1, E1max = dEDet->GetEIncOfMaxDeltaE(part.GetZ(),part.GetA());
E1 = (E1min + E1max) / 2.;
while ((E1max - E1min) > 0.1)
@@ -1515,6 +1515,7 @@
E2 = (E2max + E2min) / 2.;
}
}
+<<<<<<< TREE
Info("MakeEDeltaEZGrid","Z= %d, E2=%lf",z,E2);
@@ -1527,6 +1528,9 @@
E2 = dEDet->GetEmaxVedaloss(z)*part.GetA();
}
+=======
+
+>>>>>>> MERGE-SOURCE
KVIDentifier *line = NewLine("ID");
Add("ID", line);
line->SetZ(z);
=== modified file 'KVMultiDet/idtelescopes/KVIDTelescope.cpp'
--- KVMultiDet/idtelescopes/KVIDTelescope.cpp 2011-01-20 13:05:37 +0000
+++ KVMultiDet/idtelescopes/KVIDTelescope.cpp 2011-03-14 15:58:00 +0000
@@ -674,6 +674,7 @@
if( e1< 0.0 ) e1 = 0.0;
else{
d1->SetEnergyLoss(e1);
+ d1->SetEResAfterDetector(e2);
e1 = d1->GetCorrectedEnergy(z,a);
//status code
fCalibStatus = kCalibStatus_Calculated;
@@ -706,7 +707,9 @@
else if(d2){//2nd detector is calibrated too: get corrected energy loss
e2 = d2->GetCorrectedEnergy(z, a, -1, kFALSE);//N.B.: transmission=kFALSE because particle assumed to stop in d2
-
+ // recalculate corrected energy in first stage using info on Eres
+ d1->SetEResAfterDetector(e2);
+ e1 = d1->GetCorrectedEnergy(z, a);
}
//incident energy of particle (before 1st member of telescope)
@@ -736,6 +739,7 @@
//significantly larger, there may be a second particle.
e1 = det->GetDeltaEFromERes(z,a,einc);
if( e1< 0.0 ) e1 = 0.0;
+ det->SetEResAfterDetector(einc);
dE = det->GetCorrectedEnergy(z,a);
einc += dE;
}
@@ -763,6 +767,7 @@
//status code
fCalibStatus = kCalibStatus_Calculated;
}
+ det->SetEResAfterDetector(einc);
e1 = det->GetCorrectedEnergy(z,a,e1);
einc += e1;
}
@@ -840,8 +845,8 @@
Double_t E1, E2;
//find E1
- //go from 0.1 MeV to chio->GetBraggE(part.GetZ(),part.GetA()))
- Double_t E1min = 0.1, E1max = det_de->GetBraggE(zzz,aaa);
+ //go from 0.1 MeV to chio->GetEIncOfMaxDeltaE(part.GetZ(),part.GetA()))
+ Double_t E1min = 0.1, E1max = det_de->GetEIncOfMaxDeltaE(zzz,aaa);
E1 = (E1min + E1max) / 2.;
while ((E1max - E1min) > 0.1) {
@@ -864,8 +869,8 @@
}
//add point to Bragg line
- Double_t dE_B = det_de->GetBraggDE(zzz, aaa);
- Double_t E_B = det_de->GetBraggE(zzz, aaa);
+ Double_t dE_B = det_de->GetMaxDeltaE(zzz, aaa);
+ Double_t E_B = det_de->GetEIncOfMaxDeltaE(zzz, aaa);
Double_t Eres_B = det_de->GetERes(zzz, aaa, E_B);
B_line->SetPoint(npoi_bragg++, Eres_B, dE_B);
=== added directory 'KVMultiDet/stopping'
=== added file 'KVMultiDet/stopping/KVIonRangeTable.cpp'
--- KVMultiDet/stopping/KVIonRangeTable.cpp 1970-01-01 00:00:00 +0000
+++ KVMultiDet/stopping/KVIonRangeTable.cpp 2011-03-14 15:58:00 +0000
@@ -0,0 +1,41 @@
+//Created by KVClassFactory on Thu Feb 3 10:04:41 2011
+//Author: frankland,,,,
+
+#include "KVIonRangeTable.h"
+#include <TPluginManager.h>
+#include <TError.h>
+
+ClassImp(KVIonRangeTable)
+
+////////////////////////////////////////////////////////////////////////////////
+// BEGIN_HTML <!--
+/* -->
+<h2>KVIonRangeTable</h2>
+<h4>Abstract base class for calculation of range & energy loss of charged particles in matter</h4>
+<!-- */
+// --> END_HTML
+////////////////////////////////////////////////////////////////////////////////
+
+KVIonRangeTable::KVIonRangeTable()
+{
+ // Default constructor
+}
+
+KVIonRangeTable::~KVIonRangeTable()
+{
+ // Destructor
+}
+
+KVIonRangeTable* KVIonRangeTable::GetRangeTable(const Char_t* name)
+{
+ // Generates an instance of the KVIonRangeTable plugin class corresponding to given name.
+
+ TPluginHandler *ph;
+ //check and load plugin library
+ if (!(ph = LoadPlugin("KVIonRangeTable", name))) {
+ ::Error("KVIonRangeTable::GetRangeTable", "No plugin for KVIonRangeTable with name=%s found in .kvrootrc", name);
+ return 0;
+ }
+ KVIonRangeTable *irt = (KVIonRangeTable *) ph->ExecPlugin(0);
+ return irt;
+}
=== added file 'KVMultiDet/stopping/KVIonRangeTable.h'
--- KVMultiDet/stopping/KVIonRangeTable.h 1970-01-01 00:00:00 +0000
+++ KVMultiDet/stopping/KVIonRangeTable.h 2011-03-14 15:58:00 +0000
@@ -0,0 +1,102 @@
+//Created by KVClassFactory on Thu Feb 3 10:04:41 2011
+//Author: frankland,,,,
+
+#ifndef KVIONRANGETABLE_H
+#define KVIONRANGETABLE_H
+
+#include "KVBase.h"
+#include "KVUnits.h"
+
+class KVIonRangeTable : public KVBase {
+
+public:
+ enum SolType {
+ kEmax,
+ kEmin
+ };
+
+ KVIonRangeTable();
+ virtual ~KVIonRangeTable();
+
+ static KVIonRangeTable* GetRangeTable(const Char_t* name);
+
+ // Create and fill a list of all materials for which range tables exist.
+ // Each entry is a TNamed with the name and type (title) of the material.
+ virtual TObjArray* GetListOfMaterials() = 0;
+
+ // Return maximum incident energy (in MeV) for which range table is valid for given
+ // material and (Z,A) of incident ion
+ virtual Double_t GetEmaxValid(const Char_t* material, Int_t Z, Int_t A) = 0;
+
+ // Returns density (g/cm**3) of a material in the range tables
+ virtual Double_t GetDensity(const Char_t*) = 0;
+
+ // Set temperature (in degrees celsius) and pressure (in torr) for a given material.
+ // This usually has no effect except for gaseous materials, for which T & P determine the density (in g/cm**3).
+ virtual void SetTemperatureAndPressure(const Char_t*, Double_t temperature, Double_t pressure) = 0;
+
+ // Returns atomic number of a material in the range tables
+ virtual Double_t GetZ(const Char_t*) = 0;
+
+ // Returns atomic mass of a material in the range tables
+ virtual Double_t GetAtomicMass(const Char_t*) = 0;
+
+ // Return kTRUE if material is in range tables
+ virtual Bool_t IsMaterialKnown(const Char_t*) = 0;
+
+ // Return kTRUE if material is gaseous
+ virtual Bool_t IsMaterialGas(const Char_t*) = 0;
+
+ // Return name of material of given type or name if it is in range tables
+ virtual const Char_t* GetMaterialName(const Char_t*) = 0;
+
+ // Returns range (in mg/cm**2) of ion (Z,A) with energy E (MeV) in material.
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+ virtual Double_t GetRangeOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E,
+ Double_t Amat=0., Double_t T=-1., Double_t P=-1.) = 0;
+
+ // Returns linear range (in cm) of ion (Z,A) with energy E (MeV) in material.
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+ virtual Double_t GetLinearRangeOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E,
+ Double_t Amat=0., Double_t T=-1., Double_t P=-1.) = 0;
+
+ // Returns energy lost (in MeV) by ion (Z,A) with energy E (MeV) after thickness r (in mg/cm**2).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+ virtual Double_t GetDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E, Double_t r,
+ Double_t Amat=0., Double_t T=-1., Double_t P=-1.) = 0;
+
+ // Returns energy lost (in MeV) by ion (Z,A) with energy E (MeV) after thickness d (in cm).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+ virtual Double_t GetLinearDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E, Double_t d,
+ Double_t Amat=0., Double_t T=-1., Double_t P=-1.) = 0;
+
+ virtual Double_t GetEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E, Double_t r,
+ Double_t Amat=0., Double_t T=-1., Double_t P=-1.)=0;
+ virtual Double_t GetLinearEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E, Double_t d,
+ Double_t Amat=0., Double_t T=-1., Double_t P=-1.)=0;
+
+ virtual Double_t GetEIncFromEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t Eres, Double_t e, Double_t isoAmat = 0., Double_t T = -1., Double_t P = -1.)=0;
+ virtual Double_t GetLinearEIncFromEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t Eres, Double_t e, Double_t isoAmat = 0., Double_t T = -1., Double_t P = -1.)=0;
+
+ virtual Double_t GetEIncFromDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t DeltaE, Double_t e, enum SolType type = kEmax, Double_t isoAmat = 0., Double_t T = -1., Double_t P = -1.)=0;
+ virtual Double_t GetLinearEIncFromDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t DeltaE, Double_t e, enum SolType type = kEmax, Double_t isoAmat = 0., Double_t T = -1., Double_t P = -1.)=0;
+
+ virtual Double_t GetDeltaEFromEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t ERes, Double_t e, Double_t isoAmat = 0., Double_t T=-1., Double_t P=-1.)=0;
+ virtual Double_t GetLinearDeltaEFromEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t ERes, Double_t e, Double_t isoAmat = 0., Double_t T=-1., Double_t P=-1.)=0;
+
+ virtual Double_t GetPunchThroughEnergy(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.) =0;
+ virtual Double_t GetLinearPunchThroughEnergy(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.)=0 ;
+
+ virtual Double_t GetMaxDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.) =0;
+ virtual Double_t GetEIncOfMaxDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.) =0;
+ virtual Double_t GetLinearMaxDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.) =0;
+ virtual Double_t GetLinearEIncOfMaxDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.) =0;
+
+ ClassDef(KVIonRangeTable, 1) //Abstract base class for calculation of range & energy loss of charged particles in matter
+};
+
+#endif
=== added file 'KVMultiDet/stopping/KVUnits.cpp'
--- KVMultiDet/stopping/KVUnits.cpp 1970-01-01 00:00:00 +0000
+++ KVMultiDet/stopping/KVUnits.cpp 2011-03-14 15:58:00 +0000
@@ -0,0 +1,63 @@
+#include "KVUnits.h"
+
+NamespaceImp(KVUnits)
+
+////////////////////////////////////////////////////////////////////////////////
+// BEGIN_HTML <!--
+/* -->
+<h2>KVUnits</h2>
+<h4>Standard units of length, mass, and pressure, and their conversion factors</h4>
+<pre><span class="keyword"></span>namespace KVUnits {<br> // lengths<br>const Double_t cm = 1.0;<br>const Double_t um = 1.e-4;<br>const Double_t mm = 1.e-1;<br>const Double_t m = 1.e+2;<br> // masses<br>const Double_t g = 1.0;<br>const Double_t kg = 1.e+3;<br>const Double_t mg = 1.e-3;<br>const Double_t ug = 1.e-6;<br> // pressures<br>const Double_t torr = 1.0;<br>const Double_t atm = 760.;<br>const Double_t Pa = atm/101325.;<br>const Double_t mbar = 100.*Pa;<br>};<br><br></pre>
+This is a set of numerical constants used to define and convert units
+of length, mass, and pressure. The standard unit defined for each
+quantity is as follows:<br>
+<br>
+<div style="margin-left: 40px;">standard mass unit : grammes
+(KVUnits::g)<br>
+standard length unit : centimetres (KVUnits::cm)<br>
+standard pressure unit : torr (KVUnits::torr)<br>
+</div>
+<br>
+The value of the corresponding numerical constant for each of the
+standard units is 1. The other available units <br>
+are:<br>
+<br>
+<div style="margin-left: 40px;">mass units: microgramme (KVUnits::ug),
+milligramme (KVUnits::mg), kilogramme (KVUnits::kg)<br>
+length units: micron/micrometre (KVUnits::um), millimetre
+(KVUnits::mm), metre (KVUnits::m)<br>
+pressure units: millibar (KVUnits::mbar), pascal (KVUnits::Pa),
+atmosphere (KVUnits::atm)<br>
+</div>
+<h4>Unit conversion</h4>
+If x is a quantity expressed in terms of one of the standard units, the
+corresponding quantity in terms of a different unit is obtained by
+dividing x by the appropriate numerical constant:<br>
+<br>
+<div style="margin-left: 40px;">x = pressure in torr:
+(x/KVUnits::mbar) is pressure in millibar;<br>
+x = mass in grammes: (x/KVUnits::mg) is mass in milligrammes;<br>
+x = length in centimetres: (x/KVUnits::um) is length in
+microns/micrometres.<br>
+</div>
+<br>
+On the other hand, if x is a quantity expressed in arbitrary units,
+then in order to express it in terms of standard units multiply x by
+the appropriate numerical constant:<br>
+<br>
+<div style="margin-left: 40px;">x = pressure in atmospheres:
+(x*KVUnits::atm) is pressure in torr;<br>
+x = mass in microgrammes: (x*KVUnits::ug) is mass in grammes;<br>
+x = length in millimetres: (x*KVUnits::mm) is length in centimetres.<br>
+</div>
+<h4>Composite units</h4>
+Similar conversions can be achieved by combinations of numerical
+constants. For example, if x is a density expressed in kilogrammes per
+cubic metre, the corresponding density in standard units (g/cm**3) is<br>
+<br>
+<div style="margin-left: 40px;">x*KVUnits::kg/pow(KVUnits::m, 3)<br>
+</div><!-- */
+// --> END_HTML
+////////////////////////////////////////////////////////////////////////////////
+
+
=== added file 'KVMultiDet/stopping/KVUnits.h'
--- KVMultiDet/stopping/KVUnits.h 1970-01-01 00:00:00 +0000
+++ KVMultiDet/stopping/KVUnits.h 2011-03-14 15:58:00 +0000
@@ -0,0 +1,29 @@
+#ifndef KVUNITS__H
+#define KVUNITS__H
+
+#include <Rtypes.h>
+
+namespace KVUnits {
+ // UNITS
+ // Standard units are:
+ // [L] cm
+ // [M] g
+ // [P] Torr
+ // lengths
+const long double cm = 1.0l;
+const long double um = 1.e-4l;
+const long double mm = 1.e-1l;
+const long double m = 1.e+2l;
+ // masses
+const long double g = 1.0l;
+const long double kg = 1.e+3l;
+const long double mg = 1.e-3l;
+const long double ug = 1.e-6l;
+ // pressures
+const long double torr = 1.0l;
+const long double atm = 760.l;
+const long double Pa = atm/101325.l;
+const long double mbar = 100.l*Pa;
+};
+
+#endif
=== added file 'KVMultiDet/stopping/KVedaLoss.cpp'
--- KVMultiDet/stopping/KVedaLoss.cpp 1970-01-01 00:00:00 +0000
+++ KVMultiDet/stopping/KVedaLoss.cpp 2011-03-14 15:58:00 +0000
@@ -0,0 +1,434 @@
+//Created by KVClassFactory on Wed Feb 2 15:49:27 2011
+//Author: frankland,,,,
+
+#include "KVedaLoss.h"
+#include "KVedaLossMaterial.h"
+#include <TString.h>
+#include <TSystem.h>
+#include <TEnv.h>
+
+#define FIND_MAT_AND_EXEC(method,defval) \
+ KVedaLossMaterial* M = GetMaterial(mat); \
+ if(M) return M->method; \
+ return defval
+
+ClassImp(KVedaLoss)
+
+////////////////////////////////////////////////////////////////////////////////
+// BEGIN_HTML <!--
+/* -->
+<h2>KVedaLoss</h2>
+<h4>C++ implementation of VEDALOSS stopping power calculation</h4>
+<h3>Credits</h3>
+From the original documentation of vedaloss.f:<br>
+<pre>
+ vedaloss.doc MAJ 5/9/97
+ -------------
+ Routines ecrites par E. de Filippo
+ Mise a jour : 19/12/96 par O. LOPEZ pour UNIX sur ANASTASIE.
+ Mise a jour : 22/05/97 par J-L. Charvet (voir routine DELTA)
+ Mise a jour : 27/10/97 par R. Dayras
+ Mise a jour : 23/03/00 par R. Dayras (ajout Nobium)
+</pre>
+<!-- */
+// --> END_HTML
+////////////////////////////////////////////////////////////////////////////////
+
+KVHashList* KVedaLoss::fMaterials = 0x0;
+
+KVedaLoss::KVedaLoss()
+{
+ // Default constructor
+ SetName("VEDALOSS");
+ SetTitle("Calculation of range and energy loss of charged particles in matter using VEDALOSS range tables");
+ if (!CheckMaterialsList()) {
+ Error("KVedaLoss", "Problem reading range tables. Do not use.");
+ }
+}
+
+KVedaLoss::~KVedaLoss()
+{
+ // Destructor
+}
+
+Bool_t KVedaLoss::init_materials() const
+{
+ // PRIVATE method - called to initialize fMaterials list of all known materials
+ // properties, read from file given by TEnv variable KVedaLoss.RangeTables
+
+ Info("init_materials", "Initialising KVedaLoss...");
+ printf("\n");
+ printf("\t*************************************************************************\n");
+ printf("\t* VEDALOSS STOPPING POWER & RANGE TABLES *\n");
+ printf("\t* *\n");
+ int mat_count = 0;
+ fMaterials = new KVHashList;
+ fMaterials->SetName("VEDALOSS materials list");
+ fMaterials->SetOwner();
+
+ TString DataFilePath = gEnv->GetValue("KVedaLoss.RangeTables", "");
+ gSystem->ExpandPathName(DataFilePath);
+
+ Char_t name[25], gtype[25], state[10];
+ Float_t Amat = 0.;
+ Float_t Dens = 0.;
+ Float_t MoleWt = 0.;
+ Float_t Temp = 19.;
+ Float_t Zmat = 0.;
+
+ FILE *fp;
+ if (!(fp = fopen(DataFilePath.Data(), "r"))) {
+ Error("init_materials()", "Range tables file %s cannot be opened", DataFilePath.Data());
+ return kFALSE;
+ } else {
+ char line[132];
+ while (fgets(line, 132, fp)) { // read lines from file
+
+ switch (line[0]) {
+
+ case '/': // ignore comment lines
+ break;
+
+ case '+': // header lines
+
+ if (sscanf(line, "+ %s %s %s %f %f %f %f %f",
+ gtype, name, state, &Dens, &Zmat, &Amat,
+ &MoleWt, &Temp)
+ != 8) {
+ Error("init_materials()", "Problem reading file %s", DataFilePath.Data());
+ fclose(fp);
+ return kFALSE;
+ }
+//found a new material
+ KVedaLossMaterial *tmp_mat = new KVedaLossMaterial(name, gtype, state, Dens,
+ Zmat, Amat, MoleWt);
+ fMaterials->Add(tmp_mat);
+ if (!tmp_mat->ReadRangeTable(fp)) return kFALSE;
+ ++mat_count;
+ Double_t rho = 0.;
+ if(tmp_mat->IsGas()) tmp_mat->SetTemperatureAndPressure(19., 1.*KVUnits::atm);
+ rho = tmp_mat->GetDensity();
+ printf("\t* %2d. %-7s %-18s Z=%2d A=%5.1f rho=%6.3f g/cm**3 *\n",
+ mat_count, tmp_mat->GetType(), tmp_mat->GetName(),
+ (int)tmp_mat->GetZ(), tmp_mat->GetMass(),
+ rho);
+ break;
+ }
+ }
+ fclose(fp);
+ }
+ printf("\t* *\n");
+ printf("\t* TF1::Range::Npx = %4d TF1::EnergyLoss::Npx = %4d *\n",
+ gEnv->GetValue("KVedaLoss.Range.Npx",100), gEnv->GetValue("KVedaLoss.EnergyLoss.Npx",100));
+ printf("\t* TF1::ResidualEnergy::Npx = %4d *\n",
+ gEnv->GetValue("KVedaLoss.ResidualEnergy.Npx",100));
+ printf("\t* *\n");
+ printf("\t* INITALISATION COMPLETE *\n");
+ printf("\t*************************************************************************\n");
+ return kTRUE;
+}
+
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetAtomicMass(const Char_t* material)
+{
+ // Return atomic mass of a material in the range table.
+ // "material" can be either the type or the name of the material.
+ // Prints a warning and returns 0 if material is unknown.
+
+ KVedaLossMaterial* M = GetMaterial(material);
+ if (!M) {
+ Warning("GetAtomicMass", "Material %s is unknown. Returned mass = 0.", material);
+ return 0.0;
+ }
+ return M->GetMass();
+}
+
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetZ(const Char_t* material)
+{
+ // Return atomic number of a material in the range table.
+ // "material" can be either the type or the name of the material.
+ // Prints a warning and returns 0 if material is unknown.
+
+ KVedaLossMaterial* M = GetMaterial(material);
+ if (!M) {
+ Warning("GetZ", "Material %s is unknown. Returned Z = 0.", material);
+ return 0.0;
+ }
+ return M->GetZ();
+}
+
+//________________________________________________________________________________//
+
+Bool_t KVedaLoss::IsMaterialKnown(const Char_t* material)
+{
+ // Returns kTRUE if material of given name or type is in range table.
+ KVedaLossMaterial* M = GetMaterial(material);
+ return (M != 0x0);
+}
+
+//________________________________________________________________________________//
+
+Bool_t KVedaLoss::IsMaterialGas(const Char_t* mat)
+{
+ // Returns kTRUE if material of given name or type is gaseous.
+ FIND_MAT_AND_EXEC(IsGas(),kFALSE);
+}
+
+//________________________________________________________________________________//
+
+KVedaLossMaterial* KVedaLoss::GetMaterial(const Char_t* material)
+{
+ // Returns pointer to material of given name or type.
+ KVedaLossMaterial* M = (KVedaLossMaterial*)fMaterials->FindObject(material);
+ if (!M) {
+ M = (KVedaLossMaterial*)fMaterials->FindObjectByType(material);
+ }
+ return M;
+}
+
+//________________________________________________________________________________//
+
+const Char_t* KVedaLoss::GetMaterialName(const Char_t* mat)
+{
+ // Return name of material of given type or name if it is in range tables
+ FIND_MAT_AND_EXEC(GetName(),"");
+}
+
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetDensity(const Char_t* mat)
+{
+ // Return density of material (g/cm**3) of given type or name if it is in range tables
+ FIND_MAT_AND_EXEC(GetDensity(),0.0);
+}
+
+//________________________________________________________________________________//
+
+void KVedaLoss::SetTemperatureAndPressure(const Char_t*material, Double_t temperature, Double_t pressure)
+{
+ // Set temperature (in degrees celsius) and pressure (in torr) for a given
+ // material. This has no effect except for gaseous materials, for which T & P
+ // determine the density (in g/cm**3).
+
+ KVedaLossMaterial* M = GetMaterial(material);
+ if (M) M->SetTemperatureAndPressure(temperature, pressure);
+}
+
+//________________________________________________________________________________//
+
+void KVedaLoss::Print(Option_t*) const
+{
+ printf("KVedaLoss::%s\n%s\n", GetName(), GetTitle());
+ printf("\nEnergy loss & range tables loaded for %d materials:\n\n", fMaterials->GetEntries());
+ fMaterials->ls();
+}
+
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetRangeOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E,
+ Double_t Amat, Double_t , Double_t)
+{
+ // Returns range (in g/cm**2) of ion (Z,A) with energy E (MeV) in material.
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ FIND_MAT_AND_EXEC(GetRangeOfIon(Z, A, E, Amat),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetLinearRangeOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E,
+ Double_t Amat, Double_t T, Double_t P)
+{
+ // Returns linear range (in cm) of ion (Z,A) with energy E (MeV) in material.
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ FIND_MAT_AND_EXEC(GetLinearRangeOfIon(Z, A, E, Amat, T, P),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E, Double_t r,
+ Double_t Amat, Double_t , Double_t )
+{
+ // Returns energy lost (in MeV) by ion (Z,A) with energy E (MeV) after thickness r (in g/cm**2).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ FIND_MAT_AND_EXEC(GetDeltaEOfIon(Z, A, E, r, Amat),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetLinearDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E, Double_t d,
+ Double_t Amat, Double_t T, Double_t P)
+{
+ // Returns energy lost (in MeV) by ion (Z,A) with energy E (MeV) after thickness d (in cm).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ FIND_MAT_AND_EXEC(GetLinearDeltaEOfIon(Z, A, E, d, Amat, T, P),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E, Double_t r,
+ Double_t Amat, Double_t , Double_t )
+{
+ // Returns residual energy (in MeV) of ion (Z,A) with incident energy E (MeV) after thickness r (in g/cm**2).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ FIND_MAT_AND_EXEC(GetEResOfIon(Z, A, E, r, Amat),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetLinearEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E, Double_t d,
+ Double_t Amat, Double_t T, Double_t P)
+{
+ // Returns residual energy (in MeV) of ion (Z,A) with incident energy E (MeV) after thickness d (in cm).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ FIND_MAT_AND_EXEC(GetLinearEResOfIon(Z, A, E, d, Amat, T, P),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetEIncFromEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t Eres, Double_t e, Double_t isoAmat, Double_t , Double_t )
+{
+ // Calculates incident energy (in MeV) of an ion (Z,A) with residual energy Eres (MeV) after thickness e (in g/cm**2).
+ // Give Amat to change default (isotopic) mass of material,
+ FIND_MAT_AND_EXEC(GetEIncFromEResOfIon(Z, A, Eres, e, isoAmat),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetLinearEIncFromEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t Eres, Double_t e,
+ Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Calculates incident energy (in MeV) of an ion (Z,A) with residual energy Eres (MeV) after thickness e (in cm).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+ FIND_MAT_AND_EXEC(GetLinearEIncFromEResOfIon(Z, A, Eres, e, isoAmat, T, P),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetEIncFromDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t DeltaE, Double_t e, enum KVIonRangeTable::SolType type, Double_t isoAmat, Double_t , Double_t )
+{
+ // Calculates incident energy (in MeV) of an ion (Z,A) from energy loss DeltaE (MeV) in thickness e (in g/cm**2).
+ // Give Amat to change default (isotopic) mass of material,
+ FIND_MAT_AND_EXEC(GetEIncFromDeltaEOfIon(Z, A, DeltaE, e, type, isoAmat),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetLinearEIncFromDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t deltaE, Double_t e, enum KVIonRangeTable::SolType type,
+ Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Calculates incident energy (in MeV) of an ion (Z,A) from energy loss DeltaE (MeV) in thickness e (in cm).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+ FIND_MAT_AND_EXEC(GetLinearEIncFromDeltaEOfIon(Z, A, deltaE, e, type, isoAmat,T,P),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetDeltaEFromEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t Eres, Double_t e, Double_t isoAmat, Double_t , Double_t )
+{
+ // Calculates incident energy (in MeV) of an ion (Z,A) from energy loss DeltaE (MeV) in thickness e (in g/cm**2).
+ // Give Amat to change default (isotopic) mass of material,
+ FIND_MAT_AND_EXEC(GetDeltaEFromEResOfIon(Z, A, Eres, e, isoAmat),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetLinearDeltaEFromEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t Eres, Double_t e, Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Calculates incident energy (in MeV) of an ion (Z,A) from energy loss DeltaE (MeV) in thickness e (in cm).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+ FIND_MAT_AND_EXEC(GetLinearDeltaEFromEResOfIon(Z, A, Eres, e, isoAmat,T,P),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetPunchThroughEnergy(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat, Double_t , Double_t )
+{
+ // Calculate incident energy (in MeV) for ion (Z,A) for which the range is equal to the
+ // given thickness e (in g/cm**2). At this energy the residual energy of the ion is (just) zero,
+ // for all energies above this energy the residual energy is > 0.
+ // Give Amat to change default (isotopic) mass of material.
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+ FIND_MAT_AND_EXEC(GetPunchThroughEnergy(Z, A, e, isoAmat),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetLinearPunchThroughEnergy(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Calculate incident energy (in MeV) for ion (Z,A) for which the range is equal to the
+ // given thickness e (in cm). At this energy the residual energy of the ion is (just) zero,
+ // for all energies above this energy the residual energy is > 0.
+ // Give Amat to change default (isotopic) mass of material.
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+ FIND_MAT_AND_EXEC(GetLinearPunchThroughEnergy(Z, A, e, isoAmat,T, P),0.0);
+}
+//________________________________________________________________________________//
+
+Double_t KVedaLoss::GetMaxDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat, Double_t , Double_t )
+{
+ // Calculate maximum energy loss (in MeV) of ion (Z,A) in given thickness e (in g/cm**2).
+ // Give Amat to change default (isotopic) mass of material.
+ FIND_MAT_AND_EXEC(GetMaxDeltaEOfIon(Z, A, e, isoAmat),0.0);
+}
+
+Double_t KVedaLoss::GetEIncOfMaxDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat, Double_t , Double_t )
+{
+ // Calculate incident energy (in MeV) corresponding to maximum energy loss of ion (Z,A)
+ // in given thickness e (in g/cm**2).
+ // Give Amat to change default (isotopic) mass of material.
+
+ FIND_MAT_AND_EXEC(GetEIncOfMaxDeltaEOfIon(Z, A, e, isoAmat),0.0);
+}
+
+Double_t KVedaLoss::GetLinearMaxDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Calculate maximum energy loss (in MeV) of ion (Z,A) in given thickness e (in cm).
+ // Give Amat to change default (isotopic) mass of material.
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ FIND_MAT_AND_EXEC(GetLinearMaxDeltaEOfIon(Z, A, e, isoAmat, T, P),0.0);
+}
+
+Double_t KVedaLoss::GetLinearEIncOfMaxDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Calculate incident energy (in MeV) corresponding to maximum energy loss of ion (Z,A)
+ // in given thickness e (in cm).
+ // Give Amat to change default (isotopic) mass of material.
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ FIND_MAT_AND_EXEC(GetLinearEIncOfMaxDeltaEOfIon(Z, A, e, isoAmat, T, P),0.0);
+}
+
+Double_t KVedaLoss::GetEmaxValid(const Char_t* mat, Int_t Z, Int_t A)
+{
+ // Returns maximum energy (in MeV) for which range table is valid
+ // for given material and incident ion (Z,A)
+
+ FIND_MAT_AND_EXEC(GetEmaxValid(Z, A),0.0);
+}
+
+//______________________________________________________________________________________//
+
+TObjArray* KVedaLoss::GetListOfMaterials()
+{
+ // Create and fill a list of all materials for which range tables exist.
+ // Each entry is a TNamed with the name and type (title) of the material.
+ // User's responsibility to delete list after use (it owns its objects).
+
+ if(CheckMaterialsList()){
+ TObjArray* list = new TObjArray(fMaterials->GetEntries());
+ list->SetOwner(kTRUE);
+ TIter next(fMaterials);
+ KVedaLossMaterial* mat;
+ while( (mat = (KVedaLossMaterial*)next()) ){
+ list->Add(new TNamed(mat->GetName(), mat->GetType()));
+ }
+ return list;
+ }
+ return 0;
+}
=== added file 'KVMultiDet/stopping/KVedaLoss.h'
--- KVMultiDet/stopping/KVedaLoss.h 1970-01-01 00:00:00 +0000
+++ KVMultiDet/stopping/KVedaLoss.h 2011-03-14 15:58:00 +0000
@@ -0,0 +1,76 @@
+//Created by KVClassFactory on Wed Feb 2 15:49:27 2011
+//Author: frankland,,,,
+
+#ifndef __KVEDALOSS_H
+#define __KVEDALOSS_H
+
+#include "KVIonRangeTable.h"
+#include "KVHashList.h"
+
+class KVedaLossMaterial;
+
+class KVedaLoss : public KVIonRangeTable {
+ static KVHashList* fMaterials;// static list of all known materials
+
+ Bool_t init_materials() const;
+ Bool_t CheckMaterialsList() const {
+ if (!fMaterials) return init_materials();
+ return kTRUE;
+ };
+
+public:
+ KVedaLoss();
+ virtual ~KVedaLoss();
+
+ Bool_t IsMaterialKnown(const Char_t* material);
+ Bool_t IsMaterialGas(const Char_t* material);
+
+ KVedaLossMaterial* GetMaterial(const Char_t* material);
+ const Char_t* GetMaterialName(const Char_t* material);
+
+ virtual Double_t GetEmaxValid(const Char_t* material, Int_t Z, Int_t A);
+
+ Double_t GetDensity(const Char_t* material);
+ void SetTemperatureAndPressure(const Char_t*material, Double_t temperature, Double_t pressure);
+ Double_t GetZ(const Char_t* material);
+ Double_t GetAtomicMass(const Char_t* material);
+
+ void Print(Option_t* = "") const;
+ TObjArray* GetListOfMaterials();
+
+ virtual Double_t GetRangeOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E,
+ Double_t Amat=0., Double_t T=-1., Double_t P=-1.);
+ virtual Double_t GetLinearRangeOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E,
+ Double_t Amat=0., Double_t T=-1., Double_t P=-1.);
+
+ virtual Double_t GetDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E, Double_t r,
+ Double_t Amat=0., Double_t T=-1., Double_t P=-1.);
+ virtual Double_t GetLinearDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E, Double_t d,
+ Double_t Amat=0., Double_t T=-1., Double_t P=-1.);
+
+ virtual Double_t GetEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E, Double_t r,
+ Double_t Amat=0., Double_t T=-1., Double_t P=-1.);
+ virtual Double_t GetLinearEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t E, Double_t d,
+ Double_t Amat=0., Double_t T=-1., Double_t P=-1.);
+
+ virtual Double_t GetEIncFromEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t Eres, Double_t e, Double_t isoAmat = 0., Double_t T = -1., Double_t P = -1.);
+ virtual Double_t GetLinearEIncFromEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t Eres, Double_t e, Double_t isoAmat = 0., Double_t T = -1., Double_t P = -1.);
+
+ virtual Double_t GetEIncFromDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t DeltaE, Double_t e, enum KVIonRangeTable::SolType type = KVIonRangeTable::kEmax, Double_t isoAmat = 0., Double_t T = -1., Double_t P = -1.);
+ virtual Double_t GetLinearEIncFromDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t DeltaE, Double_t e, enum KVIonRangeTable::SolType type = KVIonRangeTable::kEmax, Double_t isoAmat = 0., Double_t T = -1., Double_t P = -1.);
+
+ virtual Double_t GetDeltaEFromEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t ERes, Double_t e, Double_t isoAmat = 0., Double_t T=-1., Double_t P=-1.);
+ virtual Double_t GetLinearDeltaEFromEResOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t ERes, Double_t e, Double_t isoAmat = 0., Double_t T=-1., Double_t P=-1.);
+
+ virtual Double_t GetPunchThroughEnergy(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.) ;
+ virtual Double_t GetLinearPunchThroughEnergy(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.) ;
+
+ virtual Double_t GetMaxDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.);
+ virtual Double_t GetEIncOfMaxDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.);
+ virtual Double_t GetLinearMaxDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.);
+ virtual Double_t GetLinearEIncOfMaxDeltaEOfIon(const Char_t* mat, Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.);
+
+ ClassDef(KVedaLoss, 1) //C++ implementation of VEDALOSS stopping power calculation
+};
+
+#endif
=== added file 'KVMultiDet/stopping/KVedaLossMaterial.cpp'
--- KVMultiDet/stopping/KVedaLossMaterial.cpp 1970-01-01 00:00:00 +0000
+++ KVMultiDet/stopping/KVedaLossMaterial.cpp 2011-03-14 15:58:00 +0000
@@ -0,0 +1,538 @@
+//Created by KVClassFactory on Wed Feb 2 16:13:08 2011
+//Author: frankland,,,,
+
+#include "KVedaLossMaterial.h"
+#include <TMath.h>
+#include "KVIonRangeTable.h"
+#include <TEnv.h>
+
+ClassImp(KVedaLossMaterial)
+
+////////////////////////////////////////////////////////////////////////////////
+// BEGIN_HTML <!--
+/* -->
+<h2>KVedaLossMaterial</h2>
+<h4>Description of material properties used by KVedaLoss range calculation</h4>
+<!-- */
+// --> END_HTML
+////////////////////////////////////////////////////////////////////////////////
+
+KVedaLossMaterial::KVedaLossMaterial()
+ : fState("unknown"),
+ fDens(0.),
+ fZmat(0),
+ fAmat(0),
+ fMoleWt(0),
+ fRange(0)
+{
+ // Default constructor
+ for (int i = 0; i < 100; i++) {
+ fEmin[i] = 0.0;
+ fEmax[i] = 500.0;
+ }
+}
+
+KVedaLossMaterial::KVedaLossMaterial(const Char_t* name, const Char_t* type, const Char_t* state,
+ Double_t density, Double_t Z, Double_t A, Double_t MoleWt)
+ : fState(state),
+ fDens(density),
+ fZmat(Z),
+ fAmat(A),
+ fMoleWt(MoleWt),
+ fRange(0)
+{
+ // create new material
+ SetName(name);
+ SetType(type);
+ for (int i = 0; i < 100; i++) {
+ fEmin[i] = 0.0;
+ fEmax[i] = 500.0;
+ }
+}
+
+KVedaLossMaterial::~KVedaLossMaterial()
+{
+ // Destructor
+}
+
+Bool_t KVedaLossMaterial::ReadRangeTable(FILE* fp)
+{
+ // Read Z- & A-dependent range parameters for material
+ //
+ // For each material we create 3 TF1 objects:
+ // KVedaLossMaterial:[type]:Range - gives range in g/cm**2 as a function of particle energy
+ // KVedaLossMaterial:[type]:EnergyLoss - gives dE as a function of particle energy
+ // KVedaLossMaterial:[type]:ResidualEnergy - gives energy after material (0 if particle stops)
+ //
+ // The TF1::fNpx parameter for these functions is defined by the environment variables
+ //
+ // KVedaLoss.Range.Npx: 20
+ // KVedaLoss.EnergyLoss.Npx: 50
+ // KVedaLoss.ResidualEnergy.Npx: 20
+ //
+
+ char line[132];
+
+ //look for energy limits to calculation validity
+ if (!fgets(line, 132, fp)) {
+ Warning("ReadRangeTable", "Problem reading energy limits in range table file for %s (%s)",
+ GetName(), GetType());
+ return kFALSE;
+ } else {
+ while (line[0] == 'Z') {
+ Int_t z1, z2;
+ Float_t e1, e2;
+ sscanf(line, "Z = %d,%d %f < E/A < %f MeV", &z1,
+ &z2, &e1, &e2);
+ char* tmp;
+ tmp = fgets(line, 132, fp);
+ for (int i = z1; i <= z2; i++) {
+ fEmin[i - 1] = e1;
+ fEmax[i - 1] = e2;
+ }
+ }
+ }
+
+ for (register int count = 0; count < ZMAX_VEDALOSS; count++) {
+
+ if (sscanf(line, "%lf %lf %lf %lf %lf %lf %lf %lf",
+ &fCoeff[count][0], &fCoeff[count][1],
+ &fCoeff[count][2], &fCoeff[count][3],
+ &fCoeff[count][4], &fCoeff[count][5],
+ &fCoeff[count][6], &fCoeff[count][7])
+ != 8) {
+ Error("ReadRangeTable", "problem reading coeffs 0-7 in range table for %s (%s)", GetName(), GetType());
+ return kFALSE;
+ }
+ if (!fgets(line, 132, fp)) {
+ Warning("ReadRangeTable", "file too short ??? %s (%s)", GetName(), GetType());
+ return kFALSE;
+ } else {
+ if (sscanf(line, "%lf %lf %lf %lf %lf %lf",
+ &fCoeff[count][8], &fCoeff[count][9],
+ &fCoeff[count][10], &fCoeff[count][11],
+ &fCoeff[count][12], &fCoeff[count][13])
+ != 6) {
+ Error("ReadRangeTable", "problem reading coeffs 8-13 in range table for %s (%s)", GetName(), GetType());
+ return kFALSE;
+ }
+ }
+ char* tmp;
+ tmp = fgets(line, 132, fp);
+ }
+
+ // get require Npx value from (user-defined) environment variables
+ Int_t my_npx = gEnv->GetValue("KVedaLoss.Range.Npx", 100);
+
+ fRange = new TF1(Form("KVedaLossMaterial:%s:Range", GetType()), this, &KVedaLossMaterial::RangeFunc,
+ 0., 1.e+03, 3, "KVedaLossMaterial", "RangeFunc");
+ fRange->SetNpx(my_npx);
+
+ my_npx = gEnv->GetValue("KVedaLoss.EnergyLoss.Npx", 100);
+ fDeltaE = new TF1(Form("KVedaLossMaterial:%s:EnergyLoss", GetType()), this, &KVedaLossMaterial::DeltaEFunc,
+ 0., 1.e+03, 4, "KVedaLossMaterial", "DeltaEFunc");
+ fDeltaE->SetNpx(my_npx);
+
+ my_npx = gEnv->GetValue("KVedaLoss.ResidualEnergy.Npx", 100);
+ fEres = new TF1(Form("KVedaLossMaterial:%s:ResidualEnergy", GetType()), this, &KVedaLossMaterial::EResFunc,
+ 0., 1.e+03, 4, "KVedaLossMaterial", "EResFunc");
+ fEres->SetNpx(my_npx);
+
+ return kTRUE;
+}
+
+void KVedaLossMaterial::ls(Option_t*) const
+{
+ printf("KVedaLossMaterial::%s Material type : %s State : %s\n", GetName(), GetType(), fState.Data());
+ printf("\tZ=%f A=%f ", fZmat, fAmat);
+ if (IsGas()) printf(" Mole Weight = %f g.", fMoleWt);
+ if (fDens > 0) printf(" Density = %f g/cm**3", fDens);
+ printf("\n\n");
+}
+
+Double_t KVedaLossMaterial::DeltaEFunc(Double_t* E, Double_t* Mypar)
+{
+ // Function parameterising the energy loss of charged particles in this material.
+ // The incident energy E[0] is given in MeV.
+ // The energy loss is calculated in MeV.
+ //
+ //Parameters:
+ // Mypar[0] = thickness of material in g/cm**2
+ // Mypar[1] = Z of charged particle
+ // Mypar[2] = A of charged particle
+ // Mypar[3] = isotope of material element (0 if material is isotopically pure)
+
+ // if range < thickness, particle stops: dE = E0
+ Double_t R0 = RangeFunc(E, &Mypar[1]);
+ if (R0 < Mypar[0]) {
+ return E[0];
+ }
+
+ // calculate energy loss - invert range function to find E corresponding to (R0 - thickness)
+ R0 -= Mypar[0];
+
+ // invert range function to get energy after absorber
+ Double_t dE = E[0] - fRange->GetX(R0);
+ return dE;
+ /*
+ // VEDALOSS inversion of R(E)
+ // range in mg/cm**2 for VEDALOSS
+ R0 /= (KVUnits::mg / pow(KVUnits::cm,2));
+ // get parameters for this Z
+ Int_t Z = (Int_t)Mypar[1];
+ Double_t A = Mypar[2];
+ Double_t *par = fCoeff[Z - 1];
+ Double_t ranx = TMath::Log(R0);
+ Double_t ranx1 = ranx - riso;
+ Double_t depsx;
+ if (ranx1 < arm)
+ depsx = (ranx1 - adn) / adm;
+ else {
+ depsx = 0.0;
+ for (register int j = 2; j < 7; j++)
+ depsx += par[j + 7] * TMath::Power(ranx1, (Double_t) (j - 1));
+ depsx += par[8];
+ }
+
+ const Double_t PERC = 0.02;
+
+ Double_t eps1 = depsx + TMath::Log(1 - PERC);
+ Double_t eps2 = depsx + TMath::Log(1 + PERC);
+ Double_t rap = TMath::Log((1 + PERC) / (1 - PERC));
+
+ Double_t rn1 = 0.0;
+ if (TMath::Exp(eps1) < 0.1)
+ rn1 = adm * eps1 + adn;
+ else {
+ for (register int j = 1; j < 7; j++)
+ rn1 += par[j + 1] * TMath::Power(eps1, (Double_t) (j - 1));
+ }
+ Double_t rn2 = 0.0;
+ if (TMath::Exp(eps2) < 0.1)
+ rn2 = adm * eps2 + adn;
+ else {
+ for (register int j = 1; j < 7; j++)
+ rn2 += par[j + 1] * TMath::Power(eps2, (Double_t) (j - 1));
+ }
+
+ Double_t epres = eps1 + (rap / (rn2 - rn1)) * (ranx1 - rn1);
+ epres = TMath::Exp(epres);
+ Double_t eres = A * epres;
+
+ // garde-fou - calculated energy after absorber > incident energy ?!!
+ //if(eres > E[0]) return 0.0;
+ return E[0] - eres;
+ */
+}
+
+Double_t KVedaLossMaterial::EResFunc(Double_t* E, Double_t* Mypar)
+{
+ // Function parameterising the residual energy of charged particles in this material.
+ // The incident energy E[0] is given in MeV.
+ // The residual energy is calculated in MeV.
+ //
+ //Parameters:
+ // Mypar[0] = thickness of material in g/cm**2
+ // Mypar[1] = Z of charged particle
+ // Mypar[2] = A of charged particle
+ // Mypar[3] = isotope of material element (0 if material is isotopically pure)
+
+ // if range < thickness, particle stops: Eres=0
+ Double_t R0 = RangeFunc(E, &Mypar[1]);
+ if (R0 < Mypar[0]) {
+ return 0.0;
+ }
+
+ // calculate energy after absorber - invert range function to find Eres corresponding to (R0 - thickness)
+ R0 -= Mypar[0];
+
+ // invert range function to get energy after absorber
+ return fRange->GetX(R0);
+}
+
+Double_t KVedaLossMaterial::RangeFunc(Double_t* E, Double_t* Mypar)
+{
+ // Function parameterising the range of charged particles in this material.
+ // The energy E[0] is given in MeV.
+ // The range is calculated in units of g/cm**2
+ //
+ //Parameters:
+ // Mypar[0] = Z of charged particle
+ // Mypar[1] = A of charged particle
+ // Mypar[2] = isotope of material element (0 if material is not isotopically pure)
+
+ Int_t Z = (Int_t)Mypar[0];
+ Double_t A = Mypar[1];
+ // get parameters for this Z
+ Double_t *par = fCoeff[Z - 1];
+
+ // set up polynomial
+ Double_t x1 = TMath::Log(0.1);
+ Double_t x2 = TMath::Log(0.2);
+ Double_t ran = 0.0;
+ for (register int j = 2; j < 7; j++)
+ ran += par[j + 1] * TMath::Power(x2, (Double_t)(j - 1));
+ ran += par[2];
+ Double_t y2 = ran;
+ ran = 0.0;
+ for (register int jj = 2; jj < 7; jj++)
+ ran += par[jj + 1] * TMath::Power(x1, (Double_t)(jj - 1));
+ ran += par[2];
+ Double_t y1 = ran;
+ Double_t adm = (y2 - y1) / (x2 - x1);
+ Double_t adn = (y1 - adm * x1);
+ // calculate energy loss
+ Double_t eps = E[0] / A; //energy in MeV/nucleon
+ Double_t dleps = TMath::Log(eps);
+ Double_t MatIsoR = 0.;
+ if (Mypar[2] > 0.0) MatIsoR = TMath::Log(Mypar[2] / fAmat);
+ Double_t riso = TMath::Log(A / par[1]) + MatIsoR;
+
+ if (eps < 0.1)
+ ran = adm * dleps + adn;
+ else {
+ ran = 0.0;
+ for (register int j = 2; j < 7; j++)
+ ran += par[j + 1] * TMath::Power(dleps, (Double_t)(j - 1));
+ ran += par[2];
+ }
+ ran += riso;
+
+ // range in g/cm**2
+ Double_t range = TMath::Exp(ran) * KVUnits::mg / pow(KVUnits::cm, 2);
+ return range;
+}
+
+TF1* KVedaLossMaterial::GetRangeFunction(Int_t Z, Int_t A, Double_t isoAmat)
+{
+ // Return function giving range (in centimetres) as a function of energy (in MeV) for
+ // charged particles (Z,A) in this material.
+ // If required, the isotopic mass of the material can be given.
+
+ fRange->SetParameters(Z, A, isoAmat);
+ fRange->SetRange(0., GetEmaxValid(Z,A));
+ return fRange;
+}
+
+TF1* KVedaLossMaterial::GetDeltaEFunction(Double_t e, Int_t Z, Int_t A, Double_t isoAmat)
+{
+ // Return function giving energy loss (in MeV) as a function of incident energy (in MeV) for
+ // charged particles (Z,A) traversing (or not) the thickness e (in g/cm**2) of this material.
+ // If required, the isotopic mass of the material can be given.
+
+ GetRangeFunction(Z, A, isoAmat);
+ fDeltaE->SetParameters(e, Z, A, isoAmat);
+ fDeltaE->SetRange(0., GetEmaxValid(Z,A));
+ return fDeltaE;
+}
+
+TF1* KVedaLossMaterial::GetEResFunction(Double_t e, Int_t Z, Int_t A, Double_t isoAmat)
+{
+ // Return function giving residual energy (in MeV) as a function of incident energy (in MeV) for
+ // charged particles (Z,A) traversing (or not) the thickness e (in g/cm**2) of this material.
+ // If required, the isotopic mass of the material can be given.
+
+ GetRangeFunction(Z, A, isoAmat);
+ fEres->SetParameters(e, Z, A, isoAmat);
+ fEres->SetRange(0., GetEmaxValid(Z,A));
+ return fEres;
+}
+
+void KVedaLossMaterial::PrintRangeTable(Int_t Z, Int_t A, Double_t isoAmat, Double_t units, Double_t T, Double_t P)
+{
+ // Print range of element (in g/cm**2) as a function of incident energy (in MeV).
+ // For solid elements, print also the linear range (in cm). To change the default units,
+ // set optional argument units (e.g. to have linear range in microns, call with units = KVUnits::um).
+ // For gaseous elements, give the temperature (in degrees) and the pressure (in torr)
+ // in order to print the range in terms of length units.
+
+ fRange->SetParameters(Z, A, isoAmat);
+ printf(" **** VEDALOSS Range Table ****\n\n");
+ ls();
+ printf(" Element: Z=%d A=%d\n\n", Z, A);
+ printf("\tENERGY (MeV)\t\tRANGE (g/cm**2)");
+ if (!IsGas() || (IsGas() && T > 0 && P > 0)) printf("\t\tLIN. RANGE");
+ SetTemperatureAndPressure(T, P);
+ printf("\n\n");
+ for (Double_t e = 0.1; e <= 1.e+4; e *= 10) {
+ printf("\t%10.5g\t\t%10.5g", e, fRange->Eval(e));
+ if (!IsGas() || (IsGas() && T > 0 && P > 0)) printf("\t\t\t%10.5g", fRange->Eval(e) / GetDensity() / units);
+ printf("\n");
+ }
+}
+
+Double_t KVedaLossMaterial::GetRangeOfIon(Int_t Z, Int_t A, Double_t E, Double_t isoAmat)
+{
+ // Returns range (in g/cm**2) of ion (Z,A) with energy E (MeV) in material.
+ // Give Amat to change default (isotopic) mass of material,
+ TF1* f = GetRangeFunction(Z, A, isoAmat);
+ return f->Eval(E);
+}
+
+Double_t KVedaLossMaterial::GetLinearRangeOfIon(Int_t Z, Int_t A, Double_t E, Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Returns range (in cm) of ion (Z,A) with energy E (MeV) in material.
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ SetTemperatureAndPressure(T, P);
+ if (fDens > 0) return GetRangeOfIon(Z, A, E, isoAmat) / GetDensity();
+ else return 0.;
+}
+
+Double_t KVedaLossMaterial::GetDeltaEOfIon(Int_t Z, Int_t A, Double_t E, Double_t e, Double_t isoAmat)
+{
+ // Returns energy lost (in MeV) by ion (Z,A) with energy E (MeV) after thickness e (in g/cm**2).
+ // Give Amat to change default (isotopic) mass of material,
+
+ TF1* f = GetDeltaEFunction(e, Z, A, isoAmat);
+ return f->Eval(E);
+}
+
+Double_t KVedaLossMaterial::GetLinearDeltaEOfIon(Int_t Z, Int_t A, Double_t E, Double_t e,
+ Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Returns energy lost (in MeV) by ion (Z,A) with energy E (MeV) after thickness e (in cm).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ SetTemperatureAndPressure(T, P);
+ e *= GetDensity();
+ return GetDeltaEOfIon(Z, A, E, e, isoAmat);
+}
+
+Double_t KVedaLossMaterial::GetEResOfIon(Int_t Z, Int_t A, Double_t E, Double_t e,
+ Double_t isoAmat)
+{
+ // Returns energy lost (in MeV) by ion (Z,A) with energy E (MeV) after thickness e (in g/cm**2).
+ // Give Amat to change default (isotopic) mass of material,
+ TF1* f = GetEResFunction(e, Z, A, isoAmat);
+ return f->Eval(E);
+}
+
+Double_t KVedaLossMaterial::GetLinearEResOfIon(Int_t Z, Int_t A, Double_t E, Double_t e,
+ Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Returns energy lost (in MeV) by ion (Z,A) with energy E (MeV) after thickness e (in cm).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ SetTemperatureAndPressure(T, P);
+ e *= GetDensity();
+ return GetEResOfIon(Z, A, E, e, isoAmat);
+}
+
+Double_t KVedaLossMaterial::GetEIncFromEResOfIon(Int_t Z, Int_t A, Double_t Eres, Double_t e, Double_t isoAmat)
+{
+ // Calculates incident energy (in MeV) of an ion (Z,A) with residual energy Eres (MeV) after thickness e (in g/cm**2).
+ // Give Amat to change default (isotopic) mass of material,
+ GetRangeFunction(Z, A, isoAmat);
+ Double_t R0 = fRange->Eval(Eres) + e;
+ return fRange->GetX(R0);
+}
+
+Double_t KVedaLossMaterial::GetLinearEIncFromEResOfIon(Int_t Z, Int_t A, Double_t Eres, Double_t e,
+ Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Calculates incident energy (in MeV) of an ion (Z,A) with residual energy Eres (MeV) after thickness e (in cm).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ SetTemperatureAndPressure(T, P);
+ e *= GetDensity();
+ return GetEIncFromEResOfIon(Z, A, Eres, e, isoAmat);
+}
+
+Double_t KVedaLossMaterial::GetEIncFromDeltaEOfIon(Int_t Z, Int_t A, Double_t DeltaE, Double_t e, enum KVIonRangeTable::SolType type, Double_t isoAmat)
+{
+ // Calculates incident energy (in MeV) of an ion (Z,A) from energy loss DeltaE (MeV) in thickness e (in g/cm**2).
+ // Give Amat to change default (isotopic) mass of material,
+ GetDeltaEFunction(e, Z, A, isoAmat);
+ Double_t e1,e2;
+ fDeltaE->GetRange(e1,e2);
+ switch(type){
+ case KVIonRangeTable::kEmin:
+ e2=GetEIncOfMaxDeltaEOfIon(Z,A,e,isoAmat);
+ break;
+ case KVIonRangeTable::kEmax:
+ e1=GetEIncOfMaxDeltaEOfIon(Z,A,e,isoAmat);
+ break;
+ }
+ return fDeltaE->GetX(DeltaE,e1,e2);
+}
+
+Double_t KVedaLossMaterial::GetLinearEIncFromDeltaEOfIon(Int_t Z, Int_t A, Double_t deltaE, Double_t e, enum KVIonRangeTable::SolType type,
+ Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Calculates incident energy (in MeV) of an ion (Z,A) from energy loss DeltaE (MeV) in thickness e (in cm).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ SetTemperatureAndPressure(T, P);
+ e *= GetDensity();
+ return GetEIncFromDeltaEOfIon(Z, A, deltaE, e, type, isoAmat);
+}
+
+Double_t KVedaLossMaterial::GetPunchThroughEnergy(Int_t Z, Int_t A, Double_t e, Double_t isoAmat)
+{
+ // Calculate incident energy (in MeV) for ion (Z,A) for which the range is equal to the
+ // given thickness e (in g/cm**2). At this energy the residual energy of the ion is (just) zero,
+ // for all energies above this energy the residual energy is > 0.
+ // Give Amat to change default (isotopic) mass of material.
+
+ return GetRangeFunction(Z,A,isoAmat)->GetX(e);
+}
+
+Double_t KVedaLossMaterial::GetLinearPunchThroughEnergy(Int_t Z, Int_t A, Double_t e, Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Calculate incident energy (in MeV) for ion (Z,A) for which the range is equal to the
+ // given thickness e (in cm). At this energy the residual energy of the ion is (just) zero,
+ // for all energies above this energy the residual energy is > 0.
+ // Give Amat to change default (isotopic) mass of material.
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ SetTemperatureAndPressure(T, P);
+ e *= GetDensity();
+ return GetPunchThroughEnergy(Z,A,e,isoAmat);
+}
+
+Double_t KVedaLossMaterial::GetMaxDeltaEOfIon(Int_t Z, Int_t A, Double_t e, Double_t isoAmat)
+{
+ // Calculate maximum energy loss (in MeV) of ion (Z,A) in given thickness e (in g/cm**2).
+ // Give Amat to change default (isotopic) mass of material.
+
+ return GetDeltaEFunction(e,Z,A,isoAmat)->GetMaximum();
+}
+
+Double_t KVedaLossMaterial::GetEIncOfMaxDeltaEOfIon(Int_t Z, Int_t A, Double_t e, Double_t isoAmat)
+{
+ // Calculate incident energy (in MeV) corresponding to maximum energy loss of ion (Z,A)
+ // in given thickness e (in g/cm**2).
+ // Give Amat to change default (isotopic) mass of material.
+
+ return GetDeltaEFunction(e,Z,A,isoAmat)->GetMaximumX();
+}
+
+Double_t KVedaLossMaterial::GetLinearMaxDeltaEOfIon(Int_t Z, Int_t A, Double_t e, Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Calculate maximum energy loss (in MeV) of ion (Z,A) in given thickness e (in cm).
+ // Give Amat to change default (isotopic) mass of material.
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ SetTemperatureAndPressure(T, P);
+ e *= GetDensity();
+ return GetMaxDeltaEOfIon(Z,A,e,isoAmat);
+}
+
+Double_t KVedaLossMaterial::GetLinearEIncOfMaxDeltaEOfIon(Int_t Z, Int_t A, Double_t e, Double_t isoAmat, Double_t T, Double_t P)
+{
+ // Calculate incident energy (in MeV) corresponding to maximum energy loss of ion (Z,A)
+ // in given thickness e (in cm).
+ // Give Amat to change default (isotopic) mass of material.
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+
+ SetTemperatureAndPressure(T, P);
+ e *= GetDensity();
+ return GetEIncOfMaxDeltaEOfIon(Z,A,e,isoAmat);
+}
+
=== added file 'KVMultiDet/stopping/KVedaLossMaterial.h'
--- KVMultiDet/stopping/KVedaLossMaterial.h 1970-01-01 00:00:00 +0000
+++ KVMultiDet/stopping/KVedaLossMaterial.h 2011-03-14 15:58:00 +0000
@@ -0,0 +1,125 @@
+//Created by KVClassFactory on Wed Feb 2 16:13:08 2011
+//Author: frankland,,,,
+
+#ifndef __KVEDALOSSMATERIAL_H
+#define __KVEDALOSSMATERIAL_H
+
+#include "KVBase.h"
+#include <TString.h>
+#include <Riostream.h>
+#include <TF1.h>
+#include "KVIonRangeTable.h"
+
+// maximum atomic number included in range tables
+#define ZMAX_VEDALOSS 100
+#define RTT 62.36367e+03 // cm^3.Torr.K^-1.mol^-1
+#define ZERO_KELVIN 273.15
+
+class KVedaLossMaterial : public KVBase {
+
+protected:
+ TString fState; // state of material = "solid", "liquid", "gas", "unknown"
+ Double_t fDens; // density of material in g/cm**3
+ Double_t fZmat; // atomic number of material
+ Double_t fAmat; // atomic mass of material
+ Double_t fMoleWt; // only used for gases - mass of one mole of gas in grammes
+ Double_t fEmax[ZMAX_VEDALOSS]; //[ZMAX_VEDALOSS] Z-dependent maximum energy/nucleon for calculation to be valid
+ Double_t fEmin[ZMAX_VEDALOSS]; //[ZMAX_VEDALOSS] Z-dependent minimum energy/nucleon for calculation to be valid
+ Double_t fCoeff[ZMAX_VEDALOSS][14]; //[ZMAX_VEDALOSS][14] parameters for range tables
+
+ TF1* fDeltaE; // function parameterising energy loss in material
+ TF1* fEres; // function parameterising residual energy after crossing material
+ TF1* fRange; // function parameterising range of charged particles in material
+
+ Double_t DeltaEFunc(Double_t*, Double_t*);
+ Double_t EResFunc(Double_t*, Double_t*);
+ Double_t RangeFunc(Double_t*, Double_t*);
+
+public:
+ KVedaLossMaterial();
+ KVedaLossMaterial(const Char_t* name, const Char_t* type, const Char_t* state,
+ Double_t density, Double_t Z, Double_t A, Double_t MoleWt = 0.0);
+ virtual ~KVedaLossMaterial();
+
+ Bool_t ReadRangeTable(FILE* fp);
+ Float_t GetEmaxValid(Int_t Z, Int_t A) const {
+ return ((Z <= ZMAX_VEDALOSS && Z > 0) ? A*fEmax[Z - 1] : 0.0);
+ };
+ Float_t GetEminValid(Int_t Z, Int_t A) const {
+ return ((Z <= ZMAX_VEDALOSS && Z > 0) ? A*fEmin[Z - 1] : 0.0);
+ };
+ Double_t GetMass() const {
+ return fAmat;
+ };
+ Double_t GetZ() const {
+ return fZmat;
+ };
+ Double_t GetDensity() const {
+ // return density of material in g/cm**3.
+ // for gaseous materials, you must call SetGasDensity(T,P) first in order
+ // to define the temperature and pressure of the gas.
+ return fDens;
+ };
+ void SetTemperatureAndPressure(Double_t T, Double_t P) {
+ // for a gaseous material, calculate density in g/cm**3 according to given
+ // conditions of temperature (T, in degrees celsius) and pressure (P, in Torr)
+ if (IsGas() && T > 0 && P > 0) fDens = (fMoleWt * P) / ((T + ZERO_KELVIN) * RTT);
+ };
+
+ Bool_t IsGas() const {
+ // returns kTRUE if material is in gaseous state
+ return (fState == "gas");
+ };
+
+ void ls(Option_t* = "") const;
+ void Print(Option_t* = "") const {
+ ls();
+ };
+
+ TF1* GetRangeFunction(Int_t Z, Int_t A, Double_t isoAmat = 0);
+ TF1* GetDeltaEFunction(Double_t e, Int_t Z, Int_t A, Double_t isoAmat = 0);
+ TF1* GetEResFunction(Double_t e, Int_t Z, Int_t A, Double_t isoAmat = 0);
+
+ void PrintRangeTable(Int_t Z, Int_t A, Double_t isoAmat = 0, Double_t units = KVUnits::cm, Double_t T = -1, Double_t P = -1);
+
+ virtual Double_t GetRangeOfIon(Int_t Z, Int_t A, Double_t E, Double_t isoAmat = 0.);
+ virtual Double_t GetLinearRangeOfIon(Int_t Z, Int_t A, Double_t E, Double_t isoAmat = 0, Double_t T = -1., Double_t P = -1.);
+
+ virtual Double_t GetDeltaEOfIon(Int_t Z, Int_t A, Double_t E, Double_t e, Double_t isoAmat = 0.);
+ virtual Double_t GetLinearDeltaEOfIon(Int_t Z, Int_t A, Double_t E, Double_t e, Double_t isoAmat = 0., Double_t T = -1., Double_t P = -1.);
+
+ virtual Double_t GetEResOfIon(Int_t Z, Int_t A, Double_t E, Double_t e, Double_t isoAmat = 0.);
+ virtual Double_t GetLinearEResOfIon(Int_t Z, Int_t A, Double_t E, Double_t e, Double_t isoAmat = 0., Double_t T = -1., Double_t P = -1.);
+
+ virtual Double_t GetEIncFromEResOfIon(Int_t Z, Int_t A, Double_t Eres, Double_t e, Double_t isoAmat = 0.);
+ virtual Double_t GetLinearEIncFromEResOfIon(Int_t Z, Int_t A, Double_t Eres, Double_t e, Double_t isoAmat = 0., Double_t T = -1., Double_t P = -1.);
+
+ virtual Double_t GetEIncFromDeltaEOfIon(Int_t Z, Int_t A, Double_t DeltaE, Double_t e, enum KVIonRangeTable::SolType type = KVIonRangeTable::kEmax, Double_t isoAmat = 0.);
+ virtual Double_t GetLinearEIncFromDeltaEOfIon(Int_t Z, Int_t A, Double_t DeltaE, Double_t e, enum KVIonRangeTable::SolType type = KVIonRangeTable::kEmax, Double_t isoAmat = 0., Double_t T = -1., Double_t P = -1.);
+
+ virtual Double_t GetDeltaEFromEResOfIon(Int_t Z, Int_t A, Double_t ERes, Double_t e, Double_t isoAmat = 0.) {
+ // Calculates energy loss (in MeV) of an ion (Z,A) from residual energy ERes (MeV) after thickness e (in mg/cm**2).
+ // Give Amat to change default (isotopic) mass of material,
+ return GetEIncFromEResOfIon(Z, A, ERes, e, isoAmat) - ERes;
+ }
+
+ virtual Double_t GetLinearDeltaEFromEResOfIon(Int_t Z, Int_t A, Double_t ERes, Double_t e,
+ Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.) {
+ // Calculates energy loss (in MeV) of an ion (Z,A) from residual energy ERes (MeV) after thickness e (in cm).
+ // Give Amat to change default (isotopic) mass of material,
+ // give temperature (degrees C) & pressure (torr) (T,P) for gaseous materials.
+ return GetLinearEIncFromEResOfIon(Z, A, ERes, e, isoAmat, T, P) - ERes;
+ }
+
+ virtual Double_t GetPunchThroughEnergy(Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0.) ;
+ virtual Double_t GetLinearPunchThroughEnergy(Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.) ;
+
+ virtual Double_t GetMaxDeltaEOfIon(Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0.);
+ virtual Double_t GetEIncOfMaxDeltaEOfIon(Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0.);
+ virtual Double_t GetLinearMaxDeltaEOfIon(Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.);
+ virtual Double_t GetLinearEIncOfMaxDeltaEOfIon(Int_t Z, Int_t A, Double_t e, Double_t isoAmat=0., Double_t T=-1., Double_t P=-1.);
+
+ ClassDef(KVedaLossMaterial, 1) //Description of material properties used by KVedaLoss range calculation
+};
+
+#endif
=== modified file 'etc/KaliVeda.rootrc'
--- etc/KaliVeda.rootrc 2011-03-14 08:44:04 +0000
+++ etc/KaliVeda.rootrc 2011-03-14 15:58:00 +0000
@@ -683,6 +683,13 @@
#
#
+# Plugins for calculation of range & energy loss of charge particles in matter.
+# Each class must derive from KVIonRangeTable abstract base class.
+Plugin.KVIonRangeTable: VEDALOSS KVedaLoss KVMultiDet "KVedaLoss()"
+
+# Ion range table used by KVMaterial. You can change this to use a different plugin defined as above.
+KVMaterial.IonRangeTable: VEDALOSS
+
# Default name for TEnv-format file containing real thicknesses of detectors
# in multidetector arrays defined for datasets.
# See KVMultiDetArray::SetDetectorThicknesses for file format.
@@ -690,11 +697,21 @@
# Special names for special datasets:
INDRA_camp5.KVMultiDetArray.DetectorThicknesses: EtalonThicknessCamp5.dat
+# KVPulseHeightDefect TF1::fNpx parameter for TF1 energy loss
+KVPulseHeightDefect.EnergyLoss.Npx: 20
# Z limit for PHD correction for energy in silicon detectors. Correction is applied for Z > ZminForPHDCorrection
KVSilicon.ZminForPHDCorrection: 10
-# File containing range tables for calculation of energy losses in KVMaterial
-KVMaterial.RangeTables: $(KVROOT)/KVFiles/data/kvloss.data
+# File containing range tables for calculation of energy losses in KVedaLoss
+KVedaLoss.RangeTables: $(KVROOT)/KVFiles/data/kvloss.data
+# TF1::fNpx parameter for TF1 range, energy loss, residual energy
+KVedaLoss.Range.Npx: 20
+KVedaLoss.EnergyLoss.Npx: 50
+KVedaLoss.ResidualEnergy.Npx: 20
+
+# KVDetector: TF1::fNpx parameter for TF1 energy loss & residual energy
+KVDetector.EnergyLoss.Npx: 20
+KVDetector.ResidualEnergy.Npx: 50
# LifeTime tables used by KVNDTManager
KVNDTManager.LifeTimeTable: KVLifeTimeTable
=== modified file 'html/ScanClasses.cpp'
--- html/ScanClasses.cpp 2009-10-30 14:59:37 +0000
+++ html/ScanClasses.cpp 2011-03-14 15:58:00 +0000
@@ -103,6 +103,7 @@
fClassTitles->Add( new KVBase("particles", "Particles & Nuclei"));
fClassTitles->Add( new KVBase("events", "Multiparticle Events"));
fClassTitles->Add( new KVBase("detectors", "Absorbers, Targets & Detectors"));
+ fClassTitles->Add( new KVBase("stopping", "Stopping Power & Range Tables for Heavy Ions"));
fClassTitles->Add( new KVBase("geometry", "Multidetector Geometry"));
fClassTitles->Add( new KVBase("indra", "INDRA Multidetector Array"));
fClassTitles->Add( new KVBase("data_management", "Data Management"));
Follow ups
