#include <RadDamageUtil.h>

Inheritance diagram for RadDamageUtil:

Collaboration diagram for RadDamageUtil:

Public Member Functions
	RadDamageUtil (const std::string &type, const std::string &name, const IInterface *parent)

virtual StatusCode	initialize () override

virtual StatusCode	finalize () override

virtual	~RadDamageUtil ()

StatusCode	initTools ()

const StatusCode	generateRamoMap (TH3F ramPotentialMap, InDetDD::PixelModuleDesign module)

const StatusCode	generateEfieldMap (TH1F &eFieldMap, InDetDD::PixelModuleDesign module)

StatusCode	generateEfieldMap (TH1F &eFieldMap, InDetDD::PixelModuleDesign module, double fluence, double biasVoltage, int layer, const std::string &TCAD_list, bool interpolate)

const StatusCode	generateDistanceTimeMap (TH2F &distanceMap_e, TH2F &distanceMap_h, TH1F &timeMap_e, TH1F &timeMap_h, TH2F &lorentzMap_e, TH2F &lorentzMap_h, TH1F &eFieldMap, InDetDD::PixelModuleDesign module)

const std::pair< double, double >	getTrappingTimes (double fluence) const

bool	saveDebugMaps ()

ServiceHandle< StoreGateSvc > &	evtStore ()
	The standard `StoreGateSvc` (event store) Returns (kind of) a pointer to the `StoreGateSvc`. More...

const ServiceHandle< StoreGateSvc > &	evtStore () const
	The standard `StoreGateSvc` (event store) Returns (kind of) a pointer to the `StoreGateSvc`. More...

const ServiceHandle< StoreGateSvc > &	detStore () const
	The standard `StoreGateSvc/DetectorStore` Returns (kind of) a pointer to the `StoreGateSvc`. More...

virtual StatusCode	sysInitialize () override
	Perform system initialization for an algorithm. More...

virtual StatusCode	sysStart () override
	Handle START transition. More...

virtual std::vector< Gaudi::DataHandle * >	inputHandles () const override
	Return this algorithm's input handles. More...

virtual std::vector< Gaudi::DataHandle * >	outputHandles () const override
	Return this algorithm's output handles. More...

Gaudi::Details::PropertyBase &	declareProperty (Gaudi::Property< T > &t)

Gaudi::Details::PropertyBase *	declareProperty (const std::string &name, SG::VarHandleKey &hndl, const std::string &doc, const SG::VarHandleKeyType &)
	Declare a new Gaudi property. More...

Gaudi::Details::PropertyBase *	declareProperty (const std::string &name, SG::VarHandleBase &hndl, const std::string &doc, const SG::VarHandleType &)
	Declare a new Gaudi property. More...

Gaudi::Details::PropertyBase *	declareProperty (const std::string &name, SG::VarHandleKeyArray &hndArr, const std::string &doc, const SG::VarHandleKeyArrayType &)

Gaudi::Details::PropertyBase *	declareProperty (const std::string &name, T &property, const std::string &doc, const SG::NotHandleType &)
	Declare a new Gaudi property. More...

Gaudi::Details::PropertyBase *	declareProperty (const std::string &name, T &property, const std::string &doc="none")
	Declare a new Gaudi property. More...

void	updateVHKA (Gaudi::Details::PropertyBase &)

MsgStream &	msg () const

MsgStream &	msg (const MSG::Level lvl) const

bool	msgLvl (const MSG::Level lvl) const

Static Public Member Functions
static const std::pair< double, double >	getMobility (double electricField, double temperature)

static double	getTanLorentzAngle (double electricField, double temperature, double bField, bool isHole)

Protected Member Functions
void	renounceArray (SG::VarHandleKeyArray &handlesArray)
	remove all handles from I/O resolution More...

std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void >	renounce (T &h)

void	extraDeps_update_handler (Gaudi::Details::PropertyBase &ExtraDeps)
	Add StoreName to extra input/output deps as needed. More...

Private Types
typedef ServiceHandle< StoreGateSvc >	StoreGateSvc_t

Private Member Functions
	RadDamageUtil ()

Gaudi::Details::PropertyBase &	declareGaudiProperty (Gaudi::Property< T > &hndl, const SG::VarHandleKeyType &)
	specialization for handling Gaudi::Property<SG::VarHandleKey> More...

Gaudi::Details::PropertyBase &	declareGaudiProperty (Gaudi::Property< T > &hndl, const SG::VarHandleKeyArrayType &)
	specialization for handling Gaudi::Property<SG::VarHandleKeyArray> More...

Gaudi::Details::PropertyBase &	declareGaudiProperty (Gaudi::Property< T > &hndl, const SG::VarHandleType &)
	specialization for handling Gaudi::Property<SG::VarHandleBase> More...

Gaudi::Details::PropertyBase &	declareGaudiProperty (Gaudi::Property< T > &t, const SG::NotHandleType &)
	specialization for handling everything that's not a Gaudi::Property<SG::VarHandleKey> or a <SG::VarHandleKeyArray> More...

Static Private Member Functions
static double	alpha (int n, int Nrep, double a)

static double	weighting3D (double x, double y, double z, int n, int m, int Nrep, double a, double b)

static double	weighting2D (double x, double z, double Lx, double sensorThickness)

Private Attributes
Gaudi::Property< int >	m_defaultRamo

Gaudi::Property< double >	m_betaElectrons

Gaudi::Property< double >	m_betaHoles

Gaudi::Property< bool >	m_saveDebugMaps

Gaudi::Property< double >	m_fieldScale {this, "fieldScale", 1.}

ToolHandle< EfieldInterpolator >	m_EfieldInterpolator

StoreGateSvc_t	m_evtStore
	Pointer to StoreGate (event store by default) More...

StoreGateSvc_t	m_detStore
	Pointer to StoreGate (detector store by default) More...

std::vector< SG::VarHandleKeyArray * >	m_vhka

bool	m_varHandleArraysDeclared

Detailed Description

Definition at line 32 of file RadDamageUtil.h.

Member Typedef Documentation

◆ StoreGateSvc_t

typedef ServiceHandle<StoreGateSvc> AthCommonDataStore< AthCommonMsg< AlgTool > >::StoreGateSvc_t

privateinherited

Definition at line 388 of file AthCommonDataStore.h.

Constructor & Destructor Documentation

◆ RadDamageUtil() [1/2]

RadDamageUtil::RadDamageUtil	(	const std::string &	type,
		const std::string &	name,
		const IInterface *	parent
	)

Definition at line 36 of file RadDamageUtil.cxx.

36 :

37 AthAlgTool(type, name, parent) {}

◆ ~RadDamageUtil()

RadDamageUtil::~RadDamageUtil ( )

virtualdefault

◆ RadDamageUtil() [2/2]

RadDamageUtil::RadDamageUtil ( )

private

Member Function Documentation

◆ alpha()

double RadDamageUtil::alpha	(	int	n,
		int	Nrep,
		double	a
	)

staticprivate

Definition at line 146 of file RadDamageUtil.cxx.

                                                      {
   return((2 * M_PI * n) / (Nrep * a));
 }

◆ declareGaudiProperty() [1/4]

Gaudi::Details::PropertyBase& AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty	(	Gaudi::Property< T > &	hndl,
		const SG::VarHandleKeyArrayType &
	)

inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleKeyArray>

Definition at line 170 of file AthCommonDataStore.h.

   {
     return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
                                                        hndl.value(), 
                                                        hndl.documentation());
  
   }

◆ declareGaudiProperty() [2/4]

Gaudi::Details::PropertyBase& AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty	(	Gaudi::Property< T > &	hndl,
		const SG::VarHandleKeyType &
	)

inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleKey>

Definition at line 156 of file AthCommonDataStore.h.

   {
     return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
                                                        hndl.value(), 
                                                        hndl.documentation());
  
   }

◆ declareGaudiProperty() [3/4]

Gaudi::Details::PropertyBase& AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty	(	Gaudi::Property< T > &	hndl,
		const SG::VarHandleType &
	)

inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleBase>

Definition at line 184 of file AthCommonDataStore.h.

   {
     return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
                                                        hndl.value(), 
                                                        hndl.documentation());
   }

◆ declareGaudiProperty() [4/4]

Gaudi::Details::PropertyBase& AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty	(	Gaudi::Property< T > &	t,
		const SG::NotHandleType &
	)

inlineprivateinherited

specialization for handling everything that's not a Gaudi::Property<SG::VarHandleKey> or a <SG::VarHandleKeyArray>

Definition at line 199 of file AthCommonDataStore.h.

   {
     return PBASE::declareProperty(t);
   }

◆ declareProperty() [1/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty	(	const std::string &	name,
		SG::VarHandleBase &	hndl,
		const std::string &	doc,
		const SG::VarHandleType &
	)

inlineinherited

Declare a new Gaudi property.

Parameters

name	Name of the property.
hndl	Object holding the property value.
doc	Documentation string for the property.

This is the version for types that derive from SG::VarHandleBase. The property value object is put on the input and output lists as appropriate; then we forward to the base class.

Definition at line 245 of file AthCommonDataStore.h.

   {
     this->declare(hndl.vhKey());
     hndl.vhKey().setOwner(this);
  
     return PBASE::declareProperty(name,hndl,doc);
   }

◆ declareProperty() [2/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty	(	const std::string &	name,
		SG::VarHandleKey &	hndl,
		const std::string &	doc,
		const SG::VarHandleKeyType &
	)

inlineinherited

Declare a new Gaudi property.

Parameters

name	Name of the property.
hndl	Object holding the property value.
doc	Documentation string for the property.

This is the version for types that derive from SG::VarHandleKey. The property value object is put on the input and output lists as appropriate; then we forward to the base class.

Definition at line 221 of file AthCommonDataStore.h.

   {
     this->declare(hndl);
     hndl.setOwner(this);
  
     return PBASE::declareProperty(name,hndl,doc);
   }

◆ declareProperty() [3/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty	(	const std::string &	name,
		SG::VarHandleKeyArray &	hndArr,
		const std::string &	doc,
		const SG::VarHandleKeyArrayType &
	)

inlineinherited

Definition at line 259 of file AthCommonDataStore.h.

   {
  
     // std::ostringstream ost;
     // ost << Algorithm::name() << " VHKA declareProp: " << name 
     //     << " size: " << hndArr.keys().size() 
     //     << " mode: " << hndArr.mode() 
     //     << "  vhka size: " << m_vhka.size()
     //     << "\n";
     // debug() << ost.str() << endmsg;
  
     hndArr.setOwner(this);
     m_vhka.push_back(&hndArr);
  
     Gaudi::Details::PropertyBase* p =  PBASE::declareProperty(name, hndArr, doc);
     if (p != 0) {
       p->declareUpdateHandler(&AthCommonDataStore<PBASE>::updateVHKA, this);
     } else {
       ATH_MSG_ERROR("unable to call declareProperty on VarHandleKeyArray " 
                     << name);
     }
  
     return p;
  
   }

◆ declareProperty() [4/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty	(	const std::string &	name,
		T &	property,
		const std::string &	doc,
		const SG::NotHandleType &
	)

inlineinherited

Declare a new Gaudi property.

Parameters

name	Name of the property.
property	Object holding the property value.
doc	Documentation string for the property.

This is the generic version, for types that do not derive from SG::VarHandleKey. It just forwards to the base class version of declareProperty.

Definition at line 333 of file AthCommonDataStore.h.

   {
     return PBASE::declareProperty(name, property, doc);
   }

◆ declareProperty() [5/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty	(	const std::string &	name,
		T &	property,
		const std::string &	doc = `"none"`
	)

inlineinherited

Declare a new Gaudi property.

Parameters

name	Name of the property.
property	Object holding the property value.
doc	Documentation string for the property.

This dispatches to either the generic declareProperty or the one for VarHandle/Key/KeyArray.

Definition at line 352 of file AthCommonDataStore.h.

   {
     typedef typename SG::HandleClassifier<T>::type htype;
     return declareProperty (name, property, doc, htype());
   }

◆ declareProperty() [6/6]

Gaudi::Details::PropertyBase& AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( Gaudi::Property< T > & t )

inlineinherited

Definition at line 145 of file AthCommonDataStore.h.

                                                                  {
     typedef typename SG::HandleClassifier<T>::type htype;
     return AthCommonDataStore<PBASE>::declareGaudiProperty(t, htype());
   }

◆ detStore()

const ServiceHandle<StoreGateSvc>& AthCommonDataStore< AthCommonMsg< AlgTool > >::detStore ( ) const

inlineinherited

The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 95 of file AthCommonDataStore.h.

95 { return m_detStore; }

◆ evtStore() [1/2]

ServiceHandle<StoreGateSvc>& AthCommonDataStore< AthCommonMsg< AlgTool > >::evtStore ( )

inlineinherited

The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 85 of file AthCommonDataStore.h.

85 { return m_evtStore; }

◆ evtStore() [2/2]

const ServiceHandle<StoreGateSvc>& AthCommonDataStore< AthCommonMsg< AlgTool > >::evtStore ( ) const

inlineinherited

The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 90 of file AthCommonDataStore.h.

90 { return m_evtStore; }

◆ extraDeps_update_handler()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::extraDeps_update_handler ( Gaudi::Details::PropertyBase & ExtraDeps )

protectedinherited

Add StoreName to extra input/output deps as needed.

use the logic of the VarHandleKey to parse the DataObjID keys supplied via the ExtraInputs and ExtraOuputs Properties to add the StoreName if it's not explicitly given

◆ finalize()

StatusCode RadDamageUtil::finalize ( )

overridevirtual

Definition at line 460 of file RadDamageUtil.cxx.

                                    {
   ATH_MSG_DEBUG("RadDamageUtil::finalize()");
   return StatusCode::SUCCESS;
 }

◆ generateDistanceTimeMap()

const StatusCode RadDamageUtil::generateDistanceTimeMap	(	TH2F *&	distanceMap_e,
		TH2F *&	distanceMap_h,
		TH1F *&	timeMap_e,
		TH1F *&	timeMap_h,
		TH2F *&	lorentzMap_e,
		TH2F *&	lorentzMap_h,
		TH1F *&	eFieldMap,
		InDetDD::PixelModuleDesign *	module
	)

Definition at line 283 of file RadDamageUtil.cxx.

                                                                                                             {
   // Implementation for precomputed maps
   //https://gitlab.cern.ch/radiationDamageDigitization/radDamage_athena_rel22/blob/rel22_radDamageDev_master/scripts/SaveMapsForAthena.C
   //TODO: From DB call each time
   double temperature = 300;
   double bField = 2*m_fieldScale;//Tesla
   //From PixelModuleDesign: TODO
   //FIXME workaround, if PixelModuleDesign not available: retrieve sensor thickness from E field
   //double sensorThickness = module->thickness() * 1000.0;//default is 200;
   double sensorThickness = 0.2; //mm
  
   //Check if x axis (sensor depth) of E field larger than IBL sensors
   if (eFieldMap->GetXaxis()->GetXmax() > 210) { //Efield in um
     sensorThickness = 0.250;
   }
   if (module) {
     sensorThickness = module->thickness() * 1000.0;//default is 200;
   }
  
   //Y-axis is time charge carrier travelled for,
   //X-axis is initial position of charge carrier,
   //Z-axis is final position of charge carrier
   distanceMap_e = new TH2F("edistance", "Electron Distance Map", 100, 0, sensorThickness, 1000, 0, 1000); //mm by ns
   distanceMap_h = new TH2F("hdistance", "Holes Distance Map", 100, 0, sensorThickness, 1000, 0, 1000);
   //Initalize distance maps
   for (int i = 1; i <= distanceMap_e->GetNbinsX(); i++) {
     for (int j = 1; j <= distanceMap_e->GetNbinsY(); j++) {
       distanceMap_h->SetBinContent(i, j, sensorThickness); //if you travel long enough, you will reach the electrode.
       distanceMap_e->SetBinContent(i, j, 0.);              //if you travel long enough, you will reach the electrode.
     }
   }
  
   //From a given place in the sensor bulk, show time-to-electrode
   timeMap_e = new TH1F("etimes", "Electron Time Map", 100, 0, sensorThickness); //mm
   timeMap_h = new TH1F("htimes", "Hole Time Map", 100, 0, sensorThickness); //mm
  
   //X axis is initial position of charge carrier (in z)
   //Y axis is distance travelled in z by charge carrier
   //Z axis is tan( lorentz_angle )
   lorentzMap_e = new TH2F("lorentz_map_e", "Lorentz Map e", 100, 0, sensorThickness, 100, 0, sensorThickness); //mm by mm
   lorentzMap_h = new TH2F("lorentz_map_h", "Lorentz Map h", 100, 0, sensorThickness, 100, 0, sensorThickness); //mm by mm
   ATH_MSG_DEBUG("Did not find time and/or distance maps.  Will compute them from the E-field map..");
  
   for (int i = 1; i <= distanceMap_e->GetNbinsX(); i++) { //Loop over initial position of charge carrier (in z)
     double time_e = 0.; //ns
     double time_h = 0.; //ns
     double distanceTravelled_e = 0; //mm
     double distanceTravelled_h = 0; //mm
     double drift_e = 0.; //mm
     double drift_h = 0.; //mm
  
     for (int j = i; j >= 1; j--) { //Lower triangle
       double dz = distanceMap_e->GetXaxis()->GetBinWidth(j); //mm
       double z_j = distanceMap_e->GetXaxis()->GetBinCenter(j); //mm
       //printf("\n \n distance map bin center i/j: %f/%f width %f (mm) \n sensor thickness: %f \n E field value: %f in
       // bin %f \n ", z_i, z_j, dz, sensorThickness, Ez, z_i*1000);
       //
       double Ez = eFieldMap->GetBinContent(eFieldMap->GetXaxis()->FindBin(z_j * 1000)) / 1e7; // in MV/mm;
       std::pair<double, double> mu = getMobility(Ez, temperature); //mm^2/MV*ns
       if (Ez > 0) {
         //Electrons
         //double tanLorentzAngle = mu.first*bField*(1.0E-3); //rad, unit conversion; pixelPitch_eta-Field is in T =
         // V*s/m^2
         double tanLorentzAngle = getTanLorentzAngle(Ez, temperature, bField, false);
         time_e += dz / (mu.first * Ez); //mm * 1/(mm/ns) = ns
  
         //Fill: time charge carrier travelled for, given staring position (i) and final position (z_j)
         distanceMap_e->SetBinContent(i, distanceMap_e->GetYaxis()->FindBin(time_e), z_j);
  
         drift_e += dz * tanLorentzAngle; //Increment the total drift parallel to plane of sensor
         distanceTravelled_e += dz; //mm (travelled in z)
         lorentzMap_e->SetBinContent(i, j, drift_e / distanceTravelled_e);
       }
       timeMap_e->SetBinContent(i, time_e);
     }
     //Mainly copied from l416 ff changed naming k=>j
     //https://gitlab.cern.ch/radiationDamageDigitization/radDamage_athena_rel22/blob/rel22_radDamageDev_master/scripts/SaveMapsForAthena.C
     for (int j = i; j <= distanceMap_e->GetNbinsX(); j++) { //holes go the opposite direction as electrons.
       double dz = distanceMap_e->GetXaxis()->GetBinWidth(j); //similar to _h
       //double Ez = eFieldMap->GetBinContent(eFieldMap->GetXaxis()->FindBin(z_i*1000))/1e7; // in MV/mm;
       double z_j = distanceMap_e->GetXaxis()->GetBinCenter(j); //mm //similar to _h
       double Ez = eFieldMap->GetBinContent(eFieldMap->GetXaxis()->FindBin(z_j * 1000)) / 1e7; // in MV/mm;
       std::pair<double, double> mu = getMobility(Ez, temperature); //mm^2/MV*ns
       if (Ez > 0) {
         //Holes
         //std::pair<double, double> mu = getMobility(Ez, temperature); //mm^2/MV*ns
         //double tanLorentzAngle = mu.second*bField*(1.0E-3); //rad
         double tanLorentzAngle = getTanLorentzAngle(Ez, temperature, bField, true);
         time_h += dz / (mu.second * Ez); //mm * 1/(mm/ns) = ns
         distanceMap_h->SetBinContent(i, distanceMap_h->GetYaxis()->FindBin(time_h), z_j);
  
         drift_h += dz * tanLorentzAngle;
         distanceTravelled_h += dz; //mm
         lorentzMap_h->SetBinContent(i, j, drift_h / distanceTravelled_h);
       }
       timeMap_h->SetBinContent(i, time_h);
     }
   }
  
   return StatusCode::SUCCESS;
   //Finally, we make maps of the average collected charge, in order to make charge chunking corrections later.
   //TODO: talk to Ben about the above comment. Where is code?
 }

◆ generateEfieldMap() [1/2]

const StatusCode RadDamageUtil::generateEfieldMap	(	TH1F *&	eFieldMap,
		InDetDD::PixelModuleDesign *	module
	)

Definition at line 214 of file RadDamageUtil.cxx.

                                                                                                     {
   //TODO: from DB
   double biasVoltage = 600.;
   double sensorThickness = module->thickness(); //default should be 0.2?
   double fluence = 8.;//*e14 neq/cm^2
  
   eFieldMap = new TH1F("hefieldz", "hefieldz", 200, 0, sensorThickness * 1e3);
  
   std::string TCAD_list = PathResolver::find_file("ibl_TCAD_EfieldProfiles.txt", "DATAPATH"); //IBL layer
   if (sensorThickness > 0.2) {
     //is blayer
     TCAD_list = PathResolver::find_file("blayer_TCAD_EfieldProfiles.txt", "DATAPATH"); //B layer
     if (sensorThickness > 0.25) {
       ATH_MSG_WARNING("Sensor thickness (" << sensorThickness << "mm) does not match geometry provided in samples");
       return StatusCode::FAILURE;
     }
   }
  
   CHECK(m_EfieldInterpolator->loadTCADlist(TCAD_list));
   eFieldMap = (TH1F*) m_EfieldInterpolator->getEfield(fluence, biasVoltage);
   return StatusCode::SUCCESS;
 }

◆ generateEfieldMap() [2/2]

StatusCode RadDamageUtil::generateEfieldMap	(	TH1F *&	eFieldMap,
		InDetDD::PixelModuleDesign *	module,
		double	fluence,
		double	biasVoltage,
		int	layer,
		const std::string &	TCAD_list,
		bool	interpolate
	)

Definition at line 237 of file RadDamageUtil.cxx.

                                                                                                                          {
   TString id;
   //TODO adapt saving location for documentation of E field interpolation
   TString predirname = "";
  
   if (interpolate) {
     id = "Interpolation";
   } else {
     id = "TCAD";
   }
   if (interpolate) {
     CHECK(m_EfieldInterpolator->loadTCADlist(TCAD_list));
     eFieldMap = (TH1F*) m_EfieldInterpolator->getEfield(fluence, biasVoltage);
   } else {
     //retrieve E field directly from file (needs to be .dat file with table)
     CHECK(m_EfieldInterpolator->loadTCADlist(TCAD_list));
     ATH_MSG_INFO("Load Efield map from " << TCAD_list);
   }
   if (!eFieldMap) {
     ATH_MSG_ERROR("E field has not been created!");
     return StatusCode::FAILURE;
   }
   //For debugging save map
   //TCAD_list.ReplaceAll(".txt","_map.root");
   TString dirname = "layer";
   dirname += layer;
   dirname += "_fl";
   dirname += TString::Format("%.1f", fluence / (float) (1e14));
   dirname += "e14_U";
   dirname += TString::Format("%.0f", biasVoltage);
   dirname += id;
   dirname.ReplaceAll(".", "-");
   dirname = predirname + dirname;
   dirname += (".root");
   eFieldMap->SaveAs(dirname.Data(), "");
   return StatusCode::SUCCESS;
 }

◆ generateRamoMap()

const StatusCode RadDamageUtil::generateRamoMap	(	TH3F *	ramPotentialMap,
		InDetDD::PixelModuleDesign *	module
	)

Definition at line 58 of file RadDamageUtil.cxx.

                                                                                                         {
   //TODO: this needs to come from DB
   double pitchX = 0.05;
   double pitchY = 0.25;
   //TODO: from PixelModuleDesign
   double sensorThickness = module->thickness() * 1000.0;//default 200;
  
   //Ramo potential evaluated up to 2x pitch away
   //from the center of the primary pixel.
   double xmin = 0.;
   double xmax = 2 * pitchX * 1000;
   double ymin = 0.;
   double ymax = 2 * pitchY * 1000;
  
   //One bin per 10 microns.
   ramoPotentialMap = new TH3F("hramomap1", "hramomap1", ((xmax - xmin) / 10.), xmin, xmax, ((ymax - ymin) / 10.), ymin,
                               ymax, int(sensorThickness * 1000) / 10, 0., sensorThickness * 1000.);
  
   //******************
   //*** Loop in z ***
   //******************
   for (int k = 1; k <= ramoPotentialMap->GetNbinsZ(); k++) {
     //use the lower bin edge.
     double z = ramoPotentialMap->GetZaxis()->GetBinCenter(k) - ramoPotentialMap->GetZaxis()->GetBinWidth(k) / 2.;
  
     //******************
     //*** Loop in x,y ***
     //******************
     for (int i = 1; i <= ramoPotentialMap->GetNbinsX(); i++) { //Loop over x
       for (int j = 1; j <= ramoPotentialMap->GetNbinsY(); j++) { //loop over y
         double x = ramoPotentialMap->GetXaxis()->GetBinCenter(i) - ramoPotentialMap->GetXaxis()->GetBinWidth(i) / 2.;
         double y = ramoPotentialMap->GetYaxis()->GetBinCenter(j) - ramoPotentialMap->GetYaxis()->GetBinWidth(j) / 2.;
  
         //*******************************
         //*** Option A: 1D approx. in z
         //*******************************
         if (m_defaultRamo == -1) {
           if (x > (pitchX * 1000.0 * 0.5) || y > (pitchY * 1000.0 * 0.5)) {
             ramoPotentialMap->SetBinContent(i, j, k, 0.01);  //outside of the primary pixel.
             //TODO what is the last bin value? Is 0.01 the min?
           } else {
             //TODO make sure all of these eta/phi values make sense for non-barrel modules too
             //The formula below parameterises the 1D Ramo potential. See ATL-COM-INDET-2018-011 for details.
             double par_a = 3 * sensorThickness / pitchY;
             double norm = exp(-par_a) + exp(-1.0);
             double val = exp(-par_a * z / (1000.0 * sensorThickness)) + exp(-z / (1000 * sensorThickness));
             val -= norm;
             val /= (2.0 - norm);
             ramoPotentialMap->SetBinContent(i, j, k, val);
           }
         }
         //*******************************
         //*** Option B: 2D approx. in xy
         //*******************************
         else if (m_defaultRamo == 0) {
           double par_a = 10.0;
           double norm = exp(-par_a) + exp(-1.);
           double val = exp(-par_a * z / (1000 * sensorThickness)) + exp(-z / (1000 * sensorThickness));
           val -= norm;
           val /= (2. - norm);
           //From equation 16 in the RadDamageDefaults support note, using solution for weighting potential in 2D
           double productSolution = weighting2D(x, z, pitchX, sensorThickness) * weighting2D(y, z, pitchY,
                                                                                             sensorThickness) * val /
                                    (weighting2D(0, z, pitchX,
                                                 sensorThickness) * weighting2D(0, z, pitchY, sensorThickness));
           ramoPotentialMap->SetBinContent(i, j, k, productSolution);
         }
         //************************************************
         //*** Option C: Poisson's eqn. w/ simple geometry
         //************************************************
         else if (m_defaultRamo > 0) {
           double fullSolution = weighting3D(x / sensorThickness, y / sensorThickness, z / sensorThickness,
                                             m_defaultRamo, m_defaultRamo, 4, pitchX / sensorThickness,
                                             pitchY / sensorThickness); //N = 4 is arbitrary; just need something bigger
                                                                        // than ~1
           ramoPotentialMap->SetBinContent(i, j, k, fullSolution);
         }//Ramo option > 0.
       }//loop over y.
     }//loop over x.
   }//loop over z.
   return StatusCode::SUCCESS;
 } //TODO: What about debugging the ramo potential?  I vaguely remember running into issues with this.

◆ getMobility()

const std::pair< double, double > RadDamageUtil::getMobility	(	double	electricField,
		double	temperature
	)

static

Definition at line 392 of file RadDamageUtil.cxx.

                                                                                                  {
   //Returns the electron/hole mobility *in the z direction*
   //Note, this already includes the Hall scattering factor!
   //These parameterizations come from C. Jacoboni et al., Solid-State Electronics 20 89. (1977) 77.  (see also
   // https://cds.cern.ch/record/684187/files/indet-2001-004.pdf).
   //electrons
   double vsat_e = 15.3 * pow(temperature, -0.87);// mm/ns
   double ecrit_e = 1.01E-7 * pow(temperature, 1.55);// MV/mm
   double beta_e = 2.57E-2 * pow(temperature, 0.66);//dimensionless
   double r_e = 1.13 + 0.0008 * (temperature - 273.);//Hall scaling factor
   //holes
   double vsat_h = 1.62 * pow(temperature, -0.52);// mm/ns
   double ecrit_h = 1.24E-7 * pow(temperature, 1.68);// MV/mm
   double beta_h = 0.46 * pow(temperature, 0.17);
   double r_h = 0.72 - 0.0005 * (temperature - 273.);
  
   double num_e = vsat_e / ecrit_e;
   double den_e = pow(1 + pow((electricField / ecrit_e), beta_e), (1 / beta_e));
   double mobility_e = r_e * num_e / den_e;
  
   double num_h = vsat_h / ecrit_h;
   double den_h = pow(1 + pow((electricField / ecrit_h), beta_h), (1 / beta_h));
   double mobility_h = r_h * num_h / den_h;
  
   return std::make_pair(mobility_e, mobility_h);
 }

◆ getTanLorentzAngle()

double RadDamageUtil::getTanLorentzAngle	(	double	electricField,
		double	temperature,
		double	bField,
		bool	isHole
	)

static

Definition at line 424 of file RadDamageUtil.cxx.

                                                                                                              {
   double hallEffect = 1.;//already in mobility//= 1.13 + 0.0008*(temperature - 273.0); //Hall Scattering Factor - taken
                          // from https://cds.cern.ch/record/684187/files/indet-2001-004.pdf
  
   if (isHole) hallEffect = 0.72 - 0.0005 * (temperature - 273.0);
   std::pair<double, double> mobility = getMobility(electricField, temperature);
   double mobility_object = mobility.first;
   if (isHole) mobility_object = mobility.second;
   double tanLorentz = hallEffect * mobility_object * bField * (1.0E-3);  //unit conversion
   return tanLorentz;
 }

◆ getTrappingTimes()

const std::pair< double, double > RadDamageUtil::getTrappingTimes ( double fluence ) const

Definition at line 439 of file RadDamageUtil.cxx.

                                                                                   {
   double trappingTimeElectrons(0.), trappingTimeHoles(0.);
  
   if (fluence != 0.0) {
     trappingTimeElectrons = 1.0 / (m_betaElectrons * fluence); //Make memberVar
     trappingTimeHoles = 1.0 / (m_betaHoles * fluence); //ns
   } else {//fluence = 0 so do not trap!
     trappingTimeElectrons = 1000; //~infinity
     trappingTimeHoles = 1000;
   }
  
   return std::make_pair(trappingTimeElectrons, trappingTimeHoles);
 }

◆ initialize()

StatusCode RadDamageUtil::initialize ( )

overridevirtual

Definition at line 44 of file RadDamageUtil.cxx.

                                      {
   ATH_CHECK(AthAlgTool::initialize());
   ATH_CHECK(m_EfieldInterpolator.retrieve());
   ATH_MSG_DEBUG("You are using RadDamageUtil for solid-state silicon detectors.");
   return StatusCode::SUCCESS;
 }

◆ initTools()

StatusCode RadDamageUtil::initTools ( )

◆ inputHandles()

virtual std::vector<Gaudi::DataHandle*> AthCommonDataStore< AthCommonMsg< AlgTool > >::inputHandles ( ) const

overridevirtualinherited

Return this algorithm's input handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ msg() [1/2]

MsgStream& AthCommonMsg< AlgTool >::msg ( ) const

inlineinherited

Definition at line 24 of file AthCommonMsg.h.

                                 {
     return this->msgStream();
   }

◆ msg() [2/2]

MsgStream& AthCommonMsg< AlgTool >::msg ( const MSG::Level lvl ) const

inlineinherited

Definition at line 27 of file AthCommonMsg.h.

                                                   {
     return this->msgStream(lvl);
   }

◆ msgLvl()

bool AthCommonMsg< AlgTool >::msgLvl ( const MSG::Level lvl ) const

inlineinherited

Definition at line 30 of file AthCommonMsg.h.

                                                {
     return this->msgLevel(lvl);
   }

◆ outputHandles()

virtual std::vector<Gaudi::DataHandle*> AthCommonDataStore< AthCommonMsg< AlgTool > >::outputHandles ( ) const

overridevirtualinherited

Return this algorithm's output handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ renounce()

std::enable_if_t<std::is_void_v<std::result_of_t<decltype(&T::renounce)(T)> > && !std::is_base_of_v<SG::VarHandleKeyArray, T> && std::is_base_of_v<Gaudi::DataHandle, T>, void> AthCommonDataStore< AthCommonMsg< AlgTool > >::renounce ( T & h )

inlineprotectedinherited

Definition at line 380 of file AthCommonDataStore.h.

   {
     h.renounce();
     PBASE::renounce (h);
   }

◆ renounceArray()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::renounceArray ( SG::VarHandleKeyArray & handlesArray )

inlineprotectedinherited

remove all handles from I/O resolution

Definition at line 364 of file AthCommonDataStore.h.

                                                           {
     handlesArray.renounce();
   }

◆ saveDebugMaps()

bool RadDamageUtil::saveDebugMaps ( )

Definition at line 453 of file RadDamageUtil.cxx.

                                   {
   return m_saveDebugMaps;
 }

◆ sysInitialize()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysInitialize ( )

overridevirtualinherited

Perform system initialization for an algorithm.

We override this to declare all the elements of handle key arrays at the end of initialization. See comments on updateVHKA.

Reimplemented in DerivationFramework::CfAthAlgTool, AthCheckedComponent< AthAlgTool >, AthCheckedComponent<::AthAlgTool >, and asg::AsgMetadataTool.

◆ sysStart()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysStart ( )

overridevirtualinherited

Handle START transition.

We override this in order to make sure that conditions handle keys can cache a pointer to the conditions container.

◆ updateVHKA()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::updateVHKA ( Gaudi::Details::PropertyBase & )

inlineinherited

Definition at line 308 of file AthCommonDataStore.h.

                                                {
     // debug() << "updateVHKA for property " << p.name() << " " << p.toString() 
     //         << "  size: " << m_vhka.size() << endmsg;
     for (auto &a : m_vhka) {
       std::vector<SG::VarHandleKey*> keys = a->keys();
       for (auto k : keys) {
         k->setOwner(this);
       }
     }
   }

◆ weighting2D()

double RadDamageUtil::weighting2D	(	double	x,
		double	z,
		double	Lx,
		double	sensorThickness
	)

staticprivate

Definition at line 195 of file RadDamageUtil.cxx.

                                                                                        {
   if (z == 0) z = 0.00001; //a pathology in the definition.M_PI
  
   //scale to binsize (inputs assumed to be in mm)
   sensorThickness *= 1000.;
   Lx *= 1000;
  
   //val is set according to equation 3 in the radDamageDefaults support note
   double val =
     (TMath::Sin(M_PI * z / sensorThickness) * TMath::SinH(0.5 * M_PI * Lx / sensorThickness) /
      (TMath::CosH(M_PI * x / sensorThickness) - TMath::Cos(M_PI * z / sensorThickness) *
       TMath::CosH(0.5 * M_PI * Lx / sensorThickness)));
   if (val > 0) return TMath::ATan(val) / M_PI;
   else return TMath::ATan(val) / M_PI + 1;
 }

◆ weighting3D()

double RadDamageUtil::weighting3D	(	double	x,
		double	y,
		double	z,
		int	n,
		int	m,
		int	Nrep,
		double	a,
		double	b
	)

staticprivate

Definition at line 155 of file RadDamageUtil.cxx.

                                                                                                           {
   //TODO: talk to ben about this comment:
   //be warned that there is numerical instability if n and m are too large!  Suggest n ~ m ~ 10.
   double potential = 0.;
  
   for (int i = -n; i <= n; i++) {
     for (int j = -m; j <= m; j++) {
       double X = 0.;
       double Y = 0.;
       double Z = 0.;
       double factor_x = 0.;
       double factor_y = 0.;
       if (i == 0 && j == 0) {
         factor_x = 1. / Nrep;
         factor_y = factor_x;
         Z = 1 - z;
       } else {
         //Equation 18 & 19 in support note
         factor_x = std::sin(i * M_PI / Nrep) / (M_PI * i);
         factor_y = std::sin(j * M_PI / Nrep) / (M_PI * j);
         //Equation 17 in support note
         double norm = std::sqrt(std::pow(alpha(i, Nrep, a), 2) + std::pow(alpha(j, Nrep, b), 2));
         Z = std::sinh(norm * (1 - z)) / sinh(norm);
       }
       //Equation 18 & 19 in support note
       X = factor_x * std::cos(alpha(i, Nrep, a) * x);
       Y = factor_y * std::cos(alpha(j, Nrep, b) * y);
  
       //Equation 20 in support note
       potential += Z * X * Y;
     }
   }
   return potential;
 }

Member Data Documentation

◆ m_betaElectrons

Gaudi::Property<double> RadDamageUtil::m_betaElectrons

private

Initial value:

{
    this, "betaElectrons", 4.5e-16, "Used in trapping time calculation"
  }

Definition at line 62 of file RadDamageUtil.h.

◆ m_betaHoles

Gaudi::Property<double> RadDamageUtil::m_betaHoles

private

Initial value:

{
    this, "betaHoles", 6.0e-16, "Used in trapping time calculation"
  }

Definition at line 67 of file RadDamageUtil.h.

◆ m_defaultRamo

Gaudi::Property<int> RadDamageUtil::m_defaultRamo

private

Initial value:

{
    this, "defaultRamo", 1, "Mapping strategy of Ramo potential"
  }

Definition at line 57 of file RadDamageUtil.h.

◆ m_detStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_detStore

privateinherited

Pointer to StoreGate (detector store by default)

Definition at line 393 of file AthCommonDataStore.h.

◆ m_EfieldInterpolator

ToolHandle<EfieldInterpolator> RadDamageUtil::m_EfieldInterpolator

private

Initial value:

{
    this, "EfieldInterpolator", "EfieldInterpolator",
    "Create an Efield for fluence and bias volatge of interest based on TCAD samples"
  }

Definition at line 83 of file RadDamageUtil.h.

◆ m_evtStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_evtStore

privateinherited

Pointer to StoreGate (event store by default)

Definition at line 390 of file AthCommonDataStore.h.

◆ m_fieldScale

Gaudi::Property<double> RadDamageUtil::m_fieldScale {this, "fieldScale", 1.}

private

Definition at line 77 of file RadDamageUtil.h.

◆ m_saveDebugMaps

Gaudi::Property<bool> RadDamageUtil::m_saveDebugMaps

private

Initial value:

{
    this, "saveDebugMaps", false, "Flag to save map"
  }

Definition at line 72 of file RadDamageUtil.h.

◆ m_varHandleArraysDeclared

bool AthCommonDataStore< AthCommonMsg< AlgTool > >::m_varHandleArraysDeclared

privateinherited

Definition at line 399 of file AthCommonDataStore.h.

◆ m_vhka

std::vector<SG::VarHandleKeyArray*> AthCommonDataStore< AthCommonMsg< AlgTool > >::m_vhka

privateinherited

Definition at line 398 of file AthCommonDataStore.h.

The documentation for this class was generated from the following files:

Public Member Functions

Static Public Member Functions

Protected Member Functions

Private Types

Private Member Functions

Static Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

◆ StoreGateSvc_t

Constructor & Destructor Documentation

◆ RadDamageUtil() [1/2]

◆ ~RadDamageUtil()

◆ RadDamageUtil() [2/2]

Member Function Documentation

◆ alpha()

◆ declareGaudiProperty() [1/4]

◆ declareGaudiProperty() [2/4]

◆ declareGaudiProperty() [3/4]

◆ declareGaudiProperty() [4/4]

◆ declareProperty() [1/6]

◆ declareProperty() [2/6]

◆ declareProperty() [3/6]

◆ declareProperty() [4/6]

◆ declareProperty() [5/6]

◆ declareProperty() [6/6]

◆ detStore()

◆ evtStore() [1/2]

◆ evtStore() [2/2]

◆ extraDeps_update_handler()

◆ finalize()

◆ generateDistanceTimeMap()

◆ generateEfieldMap() [1/2]

◆ generateEfieldMap() [2/2]

◆ generateRamoMap()

◆ getMobility()

◆ getTanLorentzAngle()

◆ getTrappingTimes()

◆ initialize()

◆ initTools()

◆ inputHandles()

◆ msg() [1/2]

◆ msg() [2/2]

◆ msgLvl()

◆ outputHandles()

◆ renounce()

◆ renounceArray()

◆ saveDebugMaps()

◆ sysInitialize()

◆ sysStart()

◆ updateVHKA()

◆ weighting2D()

◆ weighting3D()

Member Data Documentation

◆ m_betaElectrons

◆ m_betaHoles

◆ m_defaultRamo

◆ m_detStore

◆ m_EfieldInterpolator

◆ m_evtStore

◆ m_fieldScale

◆ m_saveDebugMaps

◆ m_varHandleArraysDeclared

◆ m_vhka