Commit cb323f2a authored by Morfouace's avatar Morfouace
Browse files

Adding a new scorer for process: ProcessScorer

Updating Scone, ChiNu and PISTA detectors
parent 6d7e4f28
......@@ -34,6 +34,8 @@ npsimulation
NPData/*
Projects/T40/*cxx
Inputs/EnergyLoss/*.G4table
Inputs/CrossSection/G4XS*
Inputs/CrossSection/*.C
.ls_return
Documentation/user_guide.log
*.bbl
......@@ -49,4 +51,4 @@ Projects/T40/configs/*.dat
Outputs
Outputs/*
Outputs/*/*.*
NPLib/Utility/npspectra_test
\ No newline at end of file
NPLib/Utility/npspectra_test
......@@ -4,10 +4,10 @@
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 3.5 MeV
EnergyHigh= 3.5 MeV
HalfOpenAngleMin= 0 deg
HalfOpenAngleMax= 180 deg
EnergyLow= 0.5 MeV
EnergyHigh= 0.5 MeV
HalfOpenAngleMin= 25 deg
HalfOpenAngleMax= 65 deg
x0= 0 mm
y0= 0 mm
z0= 0 mm
......
......@@ -4,8 +4,8 @@
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 1. MeV
EnergyHigh= 1. MeV
EnergyLow= 10 MeV
EnergyHigh= 10 MeV
HalfOpenAngleMin= 0 deg
HalfOpenAngleMax= 180 deg
x0= 0 mm
......
......@@ -4,13 +4,13 @@
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0.1
EnergyHigh= 5.
HalfOpenAngleMin= 0
HalfOpenAngleMax= 180
EnergyLow= 100
EnergyHigh= 100
HalfOpenAngleMin= 40
HalfOpenAngleMax= 40
x0= 0
y0= 0
z0= 0
particle= alpha
particle= 14C
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Supported particle type: proton, neutron, deuton, triton, He3 , alpha
......@@ -71,9 +71,11 @@ void TPISTAPhysics::AddDetector(double R, double Theta, double Phi){
double Height = 61.8; // mm
//double Base = 95; // mm
double Base = 78.1; // mm
double NumberOfStrips = 128;
double StripPitchHeight = Height / NumberOfStrips; // mm
double StripPitchBase = Base / NumberOfStrips; // mm
double NumberOfStripsX = 122;
double NumberOfStripsY = 97;
double StripPitchHeight = Height / NumberOfStripsY; // mm
double StripPitchBase = Base / NumberOfStripsX; // mm
// Vector U on detector face (parallel to Y strips) Y strips are along X axis
......@@ -84,7 +86,6 @@ void TPISTAPhysics::AddDetector(double R, double Theta, double Phi){
TVector3 W;
// Vector C position of detector face center
TVector3 C;
C = TVector3(R*sin(Theta)*cos(Phi),
R*sin(Theta)*sin(Phi),
Height*0.5+R*cos(Theta));
......@@ -115,11 +116,11 @@ void TPISTAPhysics::AddDetector(double R, double Theta, double Phi){
Strip_1_1 = C - (0.5*Base*U + 0.5*Height*V) + U*(StripPitchBase / 2.) + V*(StripPitchHeight / 2.);
TVector3 StripPos;
for(int i=0; i<NumberOfStrips; i++){
for(int i=0; i<NumberOfStripsX; i++){
lineX.clear();
lineY.clear();
lineZ.clear();
for(int j=0; j<NumberOfStrips; j++){
for(int j=0; j<NumberOfStripsY; j++){
StripPos = Strip_1_1 + i*U*StripPitchBase + j*V*StripPitchHeight;
lineX.push_back(StripPos.X());
lineY.push_back(StripPos.Y());
......@@ -138,30 +139,35 @@ void TPISTAPhysics::AddDetector(double R, double Theta, double Phi){
///////////////////////////////////////////////////////////////////////////
TVector3 TPISTAPhysics::GetPositionOfInteraction(const int i){
TVector3 Position = TVector3(GetStripPositionX(DetectorNumber[i], StripX[i], StripY[i]),
GetStripPositionY(DetectorNumber[i], StripX[i], StripY[i]),
GetStripPositionZ(DetectorNumber[i], StripX[i], StripY[i]));
TVector3 Position = TVector3(GetStripPositionX(DetectorNumber[i], E_StripX[i], E_StripY[i]),
GetStripPositionY(DetectorNumber[i], E_StripX[i], E_StripY[i]),
GetStripPositionZ(DetectorNumber[i], E_StripX[i], E_StripY[i]));
/*TVector3 Position = TVector3(GetStripPositionX(DetectorNumber[i], DE_StripX[i], E_StripY[i]),
GetStripPositionY(DetectorNumber[i], DE_StripX[i], E_StripY[i]),
GetStripPositionZ(DetectorNumber[i], DE_StripX[i], E_StripY[i]));
*/
return Position;
}
///////////////////////////////////////////////////////////////////////////
TVector3 TPISTAPhysics::GetDetectorNormal(const int i){
TVector3 U = TVector3(GetStripPositionX(DetectorNumber[i],128,1),
GetStripPositionY(DetectorNumber[i],128,1),
GetStripPositionZ(DetectorNumber[i],128,1))
TVector3 U = TVector3(GetStripPositionX(DetectorNumber[i],122,1),
GetStripPositionY(DetectorNumber[i],122,1),
GetStripPositionZ(DetectorNumber[i],122,1))
-TVector3(GetStripPositionX(DetectorNumber[i],1,1),
GetStripPositionY(DetectorNumber[i],1,1),
GetStripPositionZ(DetectorNumber[i],1,1));
-TVector3(GetStripPositionX(DetectorNumber[i],122,1),
GetStripPositionY(DetectorNumber[i],122,1),
GetStripPositionZ(DetectorNumber[i],122,1));
TVector3 V = TVector3(GetStripPositionX(DetectorNumber[i],128,128),
GetStripPositionY(DetectorNumber[i],128,128),
GetStripPositionZ(DetectorNumber[i],128,128))
TVector3 V = TVector3(GetStripPositionX(DetectorNumber[i],122,97),
GetStripPositionY(DetectorNumber[i],122,97),
GetStripPositionZ(DetectorNumber[i],122,97))
-TVector3(GetStripPositionX(DetectorNumber[i],128,1),
GetStripPositionY(DetectorNumber[i],128,1),
GetStripPositionZ(DetectorNumber[i],128,1));
-TVector3(GetStripPositionX(DetectorNumber[i],122,1),
GetStripPositionY(DetectorNumber[i],122,1),
GetStripPositionZ(DetectorNumber[i],122,1));
TVector3 Normal = U.Cross(V);
......@@ -180,36 +186,45 @@ void TPISTAPhysics::BuildPhysicalEvent() {
PreTreat();
if(1 /*CheckEvent() == 1*/){
vector<TVector2> couple = Match_X_Y();
EventMultiplicity = couple.size();
for(unsigned int i=0; i<couple.size(); i++){
int N = m_PreTreatedData->GetFirstStage_XE_DetectorNbr(couple[i].X());
int X = m_PreTreatedData->GetFirstStage_XE_StripNbr(couple[i].X());
int Y = m_PreTreatedData->GetFirstStage_YE_StripNbr(couple[i].Y());
double XE = m_PreTreatedData->GetFirstStage_XE_Energy(couple[i].X());
double YE = m_PreTreatedData->GetFirstStage_YE_Energy(couple[i].Y());
DetectorNumber.push_back(N);
StripX.push_back(X);
StripY.push_back(Y);
DE.push_back(XE);
PosX.push_back(GetPositionOfInteraction(i).x());
PosY.push_back(GetPositionOfInteraction(i).y());
PosZ.push_back(GetPositionOfInteraction(i).z());
int SecondStageMult = m_PreTreatedData->GetSecondStageMultXEnergy();
//vector<TVector2> couple = Match_X_Y();
//EventMultiplicity = couple.size();
int FirstStageMult = m_PreTreatedData->GetFirstStageMultXEnergy();
int SecondStageMult = m_PreTreatedData->GetSecondStageMultXEnergy();
for(unsigned int i=0; i<FirstStageMult; i++){
for(unsigned int j=0; j<SecondStageMult; j++){
if(m_PreTreatedData->GetSecondStage_XE_DetectorNbr(j)==N){
double XDE = m_PreTreatedData->GetSecondStage_XE_Energy(j);
double YDE = m_PreTreatedData->GetSecondStage_YE_Energy(j);
int DE_N = m_PreTreatedData->GetFirstStage_XE_DetectorNbr(i);
int E_N = m_PreTreatedData->GetSecondStage_XE_DetectorNbr(j);
if(DE_N==E_N){
int DE_X = m_PreTreatedData->GetFirstStage_XE_StripNbr(i);
int DE_Y = m_PreTreatedData->GetFirstStage_YE_StripNbr(i);
double DE_XE = m_PreTreatedData->GetFirstStage_XE_Energy(i);
double DE_YE = m_PreTreatedData->GetFirstStage_YE_Energy(i);
int E_X = m_PreTreatedData->GetSecondStage_XE_StripNbr(j);
int E_Y = m_PreTreatedData->GetSecondStage_YE_StripNbr(j);
double E_XE = m_PreTreatedData->GetSecondStage_XE_Energy(j);
double E_YE = m_PreTreatedData->GetSecondStage_YE_Energy(j);
DetectorNumber.push_back(DE_N);
DE_StripX.push_back(DE_X);
DE_StripY.push_back(DE_Y);
DE.push_back(DE_YE);
E_StripX.push_back(E_X);
E_StripY.push_back(E_Y);
E.push_back(E_XE);
PosX.push_back(GetPositionOfInteraction(i).x());
PosY.push_back(GetPositionOfInteraction(i).y());
PosZ.push_back(GetPositionOfInteraction(i).z());
E.push_back(XDE);
}
}
}
}
EventMultiplicity = DetectorNumber.size();
}
///////////////////////////////////////////////////////////////////////////
......@@ -284,7 +299,7 @@ void TPISTAPhysics::PreTreat() {
}
}
}
unsigned int sizeBack = m_EventData->GetFirstStageMultXEnergy();
unsigned int sizeBack = m_EventData->GetFirstStageMultYEnergy();
for (UShort_t i = 0; i < sizeBack ; ++i) {
if (m_EventData->GetFirstStage_YE_Energy(i) > m_E_RAW_Threshold) {
Double_t Energy = m_EventData->GetFirstStage_YE_Energy(i);
......@@ -313,7 +328,7 @@ void TPISTAPhysics::PreTreat() {
}
}
}
sizeBack = m_EventData->GetSecondStageMultXEnergy();
sizeBack = m_EventData->GetSecondStageMultYEnergy();
for (UShort_t i = 0; i < sizeBack ; ++i) {
if (m_EventData->GetSecondStage_YE_Energy(i) > m_E_RAW_Threshold) {
Double_t Energy = m_EventData->GetSecondStage_YE_Energy(i);
......@@ -404,8 +419,10 @@ void TPISTAPhysics::Clear() {
// DSSD
DetectorNumber.clear();
E.clear();
StripX.clear();
StripY.clear();
DE_StripX.clear();
DE_StripY.clear();
E_StripX.clear();
E_StripY.clear();
Time.clear();
DE.clear();
}
......
......@@ -69,8 +69,10 @@ class TPISTAPhysics : public TObject, public NPL::VDetector {
vector<int> DetectorNumber;
vector<double> E;
vector<double> DE;
vector<int> StripX;
vector<int> StripY;
vector<int> DE_StripX;
vector<int> DE_StripY;
vector<int> E_StripY;
vector<int> E_StripX;
vector<double> Time;
vector<double> PosX;
......
......@@ -52,6 +52,13 @@ void TSconeData::Clear() {
fScone_T_DetectorNbr.clear();
fScone_T_PlasticNbr.clear();
fScone_Time.clear();
// Flah for simulation
fScone_HasCaptured.clear();
fScone_CaptureTime.clear();
fScone_GammaEnergy.clear();
fScone_ProtonEnergy.clear();
fScone_ProtonTime.clear();
fScone_FCProcess.clear();
}
......
......@@ -44,27 +44,33 @@ class TSconeData : public TObject {
vector<UShort_t> fScone_T_PlasticNbr;
vector<Double_t> fScone_Time;
//////////////////////////////////////////////////////////////
// Constructor and destructor
// Flag for simulation
vector<UShort_t> fScone_HasCaptured; //=1 for neutron capture on H, 2 on Gd 0 otherwise
vector<double> fScone_CaptureTime;
vector<double> fScone_GammaEnergy;
vector<double> fScone_ProtonEnergy;
vector<double> fScone_ProtonTime;
vector<int> fScone_FCProcess;
//////////////////////////////////////////////////////////////
// Constructor and destructor
public:
TSconeData();
~TSconeData();
//////////////////////////////////////////////////////////////
// Inherited from TObject and overriden to avoid warnings
//////////////////////////////////////////////////////////////
// Inherited from TObject and overriden to avoid warnings
public:
void Clear();
void Clear(const Option_t*) {};
void Dump() const;
//////////////////////////////////////////////////////////////
// Getters and Setters
// Prefer inline declaration to avoid unnecessary called of
// frequently used methods
// add //! to avoid ROOT creating dictionnary for the methods
//////////////////////////////////////////////////////////////
// Getters and Setters
// Prefer inline declaration to avoid unnecessary called of
// frequently used methods
// add //! to avoid ROOT creating dictionnary for the methods
public:
////////////////////// SETTERS ////////////////////////
// Energy
......@@ -82,31 +88,73 @@ class TSconeData : public TObject {
};//!
// Flag for simulation
inline void SetCapture(const UShort_t& capture){
fScone_HasCaptured.push_back(capture);
};//
inline void SetCaptureTime(double capture_time){
fScone_CaptureTime.push_back(capture_time);
};//
inline void SetGammaEnergy(double energy){
fScone_GammaEnergy.push_back(energy);
};//
inline void SetProtonEnergy(double energy){
fScone_ProtonEnergy.push_back(energy);
};//
inline void SetProtonTime(double time){
fScone_ProtonTime.push_back(time);
};//
inline void SetFCProcess(int val){
fScone_FCProcess.push_back(val);
};//
////////////////////// GETTERS ////////////////////////
// Energy
inline UShort_t GetMultEnergy() const
{return fScone_E_DetectorNbr.size();}
{return fScone_E_DetectorNbr.size();}
inline UShort_t GetE_DetectorNbr(const unsigned int &i) const
{return fScone_E_DetectorNbr[i];}//!
{return fScone_E_DetectorNbr[i];}//!
inline UShort_t GetE_PlasticNbr(const unsigned int &i) const
{return fScone_E_PlasticNbr[i];}//!
{return fScone_E_PlasticNbr[i];}//!
inline Double_t Get_Energy(const unsigned int &i) const
{return fScone_Energy[i];}//!
{return fScone_Energy[i];}//!
// Time
inline UShort_t GetMultTime() const
{return fScone_T_DetectorNbr.size();}
{return fScone_T_DetectorNbr.size();}
inline UShort_t GetT_DetectorNbr(const unsigned int &i) const
{return fScone_T_DetectorNbr[i];}//!
{return fScone_T_DetectorNbr[i];}//!
inline UShort_t GetT_PlasticNbr(const unsigned int &i) const
{return fScone_T_PlasticNbr[i];}//!
{return fScone_T_PlasticNbr[i];}//!
inline Double_t Get_Time(const unsigned int &i) const
{return fScone_Time[i];}//!
//////////////////////////////////////////////////////////////
// Required for ROOT dictionnary
ClassDef(TSconeData,1) // SconeData structure
{return fScone_Time[i];}//!
// Flag for simulation
inline UShort_t GetCapture(const unsigned int &i) const
{return fScone_HasCaptured[i];}//!
inline double GetCaptureTime(const unsigned int &i) const
{return fScone_CaptureTime[i];}//!
inline UShort_t GetGammaMult() const
{return fScone_GammaEnergy.size();}
inline double GetGammaEnergy(const unsigned int &i) const
{return fScone_GammaEnergy[i];}//!
inline UShort_t GetProtonMult() const
{return fScone_ProtonEnergy.size();}
inline double GetProtonEnergy(const unsigned int &i) const
{return fScone_ProtonEnergy[i];}//!
inline double GetProtonTime(const unsigned int &i) const
{return fScone_ProtonTime[i];}//!
inline UShort_t GetFCProcessMult() const
{return fScone_FCProcess.size();}//!
inline int GetFCProcess(const unsigned int &i) const
{return fScone_FCProcess[i];}//!
//////////////////////////////////////////////////////////////
// Required for ROOT dictionnary
ClassDef(TSconeData,1) // SconeData structure
};
#endif
......@@ -87,13 +87,24 @@ void TSconePhysics::BuildPhysicalEvent() {
// match energy and time together
unsigned int mysizeE = m_PreTreatedData->GetMultEnergy();
unsigned int mysizeT = m_PreTreatedData->GetMultTime();
vector<double> tmp_energy;
vector<int> tmp_det;
vector<int> tmp_plastic;
tmp_energy.clear();
tmp_det.clear();
tmp_plastic.clear();
//for (UShort_t e = 0; e < mysizeE ; e++) {
if (mysizeE == mysizeT) {
for (UShort_t t = 0; t < mysizeT ; t++) {
if (m_PreTreatedData->GetE_DetectorNbr(t) == m_PreTreatedData->GetT_DetectorNbr(t)) {
int det = m_PreTreatedData->GetE_DetectorNbr(t);
tmp_energy.push_back(m_PreTreatedData->Get_Energy(t));
tmp_det.push_back(m_PreTreatedData->GetE_DetectorNbr(t));
tmp_plastic.push_back(m_PreTreatedData->GetE_PlasticNbr(t));
Energy_det[det-1] += m_PreTreatedData->Get_Energy(t);
//cout << det << " " << Energy_det[det-1] << " " << m_PreTreatedData->GetE_PlasticNbr(t) << endl;
if(m_PreTreatedData->Get_Time(t)>Time_det[det-1]){
Time_det[det-1] = m_PreTreatedData->Get_Time(t);
}
......@@ -111,8 +122,33 @@ void TSconePhysics::BuildPhysicalEvent() {
DetectorNumber.push_back(i+1);
Energy.push_back(Energy_det[i]);
Time.push_back(Time_det[i]);
/*if(Energy_det[i]>15){
cout << "**********************************" << endl;
cout << "event with E<15 MeV !!!!" << endl;
double Esum=0;
for(unsigned int k=0; k<tmp_energy.size(); k++){
if(tmp_det[k]==i+1) {
cout << i+1 << " / " <<tmp_det[k] << " / " << tmp_plastic[k] << " / " << tmp_energy[k] << endl;
Esum += tmp_energy[k];
}
}
cout << "final-> " << Energy_det[i] << " / " << Esum << endl;
}*/
}
}
// Capture Flag
HasCaptured = m_PreTreatedData->GetCapture(0);
for(UShort_t i=0; i< m_PreTreatedData->GetGammaMult(); i++){
GammaEnergy.push_back(m_PreTreatedData->GetGammaEnergy(i));
}
for(UShort_t i=0; i< m_PreTreatedData->GetProtonMult(); i++){
ProtonEnergy.push_back(m_PreTreatedData->GetProtonEnergy(i));
ProtonTime.push_back(m_PreTreatedData->GetProtonTime(i));
}
}
///////////////////////////////////////////////////////////////////////////
......@@ -143,6 +179,17 @@ void TSconePhysics::PreTreat() {
Double_t Time= Cal->ApplyCalibration("Scone/TIME"+NPL::itoa(m_EventData->GetT_DetectorNbr(i)),m_EventData->Get_Time(i));
m_PreTreatedData->SetTime(m_EventData->GetT_DetectorNbr(i), m_EventData->GetT_PlasticNbr(i), Time);
}
// Capture Flag;
m_PreTreatedData->SetCapture(m_EventData->GetCapture(0));
for(UShort_t i=0; i< m_EventData->GetGammaMult(); i++){
m_PreTreatedData->SetGammaEnergy(m_EventData->GetGammaEnergy(i));
}
for(UShort_t i=0; i< m_EventData->GetProtonMult(); i++){
m_PreTreatedData->SetProtonEnergy(m_EventData->GetProtonEnergy(i));
m_PreTreatedData->SetProtonTime(m_EventData->GetProtonTime(i));
}
}
......@@ -216,6 +263,10 @@ void TSconePhysics::Clear() {
DetectorNumber.clear();
Energy.clear();
Time.clear();
HasCaptured = 0;
GammaEnergy.clear();
ProtonEnergy.clear();
ProtonTime.clear();
}
......
......@@ -65,6 +65,10 @@ class TSconePhysics : public TObject, public NPL::VDetector {
vector<int> DetectorNumber;
vector<double> Energy;
vector<double> Time;
int HasCaptured;
vector<double> GammaEnergy;
vector<double> ProtonEnergy;
vector<double> ProtonTime;
/// A usefull method to bundle all operation to add a detector
void AddDetector(TVector3 POS);
......
......@@ -288,6 +288,7 @@ void Nucleus::Extract(string line){
// life time
string s_lt_units = line.substr(69,3);
string s_LifeTime = line.substr(57,12);
string s_LifeTimeErr = line.substr(72,7);
// Remove star
replace (s_LifeTime.begin(), s_LifeTime.end(), '*' , ' ');
// Remove <
......@@ -301,45 +302,81 @@ void Nucleus::Extract(string line){
// Remove space
s_LifeTime.erase( std::remove_if( s_LifeTime.begin(), s_LifeTime.end(), ::isspace ), s_LifeTime.end() );
s_LifeTimeErr.erase( std::remove_if( s_LifeTimeErr.begin(), s_LifeTimeErr.end(), ::isspace ), s_LifeTimeErr.end() );
s_lt_units.erase( std::remove_if( s_lt_units.begin(), s_lt_units.end(), ::isspace ), s_lt_units.end() );
if(s_LifeTime=="stbl")
if(s_LifeTime=="stbl"){
fLifeTime=-1;
else
fLifeTimeErr=-1;
}
else{
fLifeTime=atof(s_LifeTime.data());
if(s_lt_units=="ms")
fLifeTimeErr=atof(s_LifeTimeErr.data());
}
if(s_lt_units=="ms"){
fLifeTime*=1e-3;
else if(s_lt_units=="us")
fLifeTimeErr*=1e-3;
}
else if(s_lt_units=="us"){
fLifeTime*=1e-6;
else if(s_lt_units=="ns")
fLifeTimeErr*=1e-6;
}
else if(s_lt_units=="ns"){
fLifeTime*=1e-9;
else if(s_lt_units=="ps")
fLifeTimeErr*=1e-9;
}
else if(s_lt_units=="ps"){
fLifeTime*=1e-12;
else if(s_lt_units=="fs")
fLifeTimeErr*=1e-12;
}
else if(s_lt_units=="fs"){
fLifeTime*=1e-15;
else if(s_lt_units=="as")
fLifeTimeErr*=1e-15;
}
else if(s_lt_units=="as"){
fLifeTime*=1e-18;
else if(s_lt_units=="zs")
fLifeTimeErr*=1e-18;
}
else if(s_lt_units=="zs"){
fLifeTime*=1e-21;
else if(s_lt_units=="ys")
fLifeTimeErr*=1e-21;
}
else if(s_lt_units=="ys"){
fLifeTime*=1e-24;
else if(s_lt_units=="m")
fLifeTimeErr*=1e-24;
}
else if(s_lt_units=="m"){
fLifeTime*=60;
else if(s_lt_units=="h")
fLifeTimeErr*=60;
}
else if(s_lt_units=="h"){
fLifeTime*=3600;
else if(s_lt_units=="d")
fLifeTimeErr*=3600;
}
else if(s_lt_units=="d"){
fLifeTime*=3600*24;
else if(s_lt_units=="y")
fLifeTimeErr*=3600*24;
}
else if(s_lt_units=="y"){
fLifeTime*=3600*24*365.25;
else if(s_lt_units=="ky")
fLifeTimeErr*=3600*24*365.25;
}
else if(s_lt_units=="ky"){
fLifeTime*=3600*24*365.25*1e3;