LArG4GenShowerLib Class Reference

Class for shower library generation algorithm. More...

#include <LArG4GenShowerLib.h>

Inheritance diagram for LArG4GenShowerLib:

Collaboration diagram for LArG4GenShowerLib:

Public Member Functions
	LArG4GenShowerLib (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	finalize ()
StatusCode	execute ()
virtual StatusCode	sysInitialize () override
	Override sysInitialize. More...
virtual const DataObjIDColl &	extraOutputDeps () const override
	Return the list of extra output dependencies. More...
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	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

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 std::map< std::string, ShowerLib::IShowerLib * >	libMap
typedef ServiceHandle< StoreGateSvc >	StoreGateSvc_t

Private Member Functions
void	clusterize (ShowerLib::StepInfoList *stepinfo)
const ShowerLib::StepInfoCollection *	getStepInfo ()
ShowerLib::StepInfoList *	copyStepInfo (const ShowerLib::StepInfoCollection *stepinfo)
ShowerLib::StepInfoList *	copyStepInfoZeroCleanup (const ShowerLib::StepInfoCollection *stepinfo)
void	truncate (ShowerLib::StepInfoList *stepinfo)
HepMC::ConstGenParticlePtr	getParticleFromMC ()
	return first MC truth particle for event More...
void	calculateMoments (const ShowerLib::StepInfoCollection &eventSteps, double &weights, double &xavfra, double &yavfra, double &ravfra)
	calculate moments from StepInfoCollection More...
void	addingTagsToLibrary ()
	adding tag information (release, detector description, ...) to library comment More...
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...

Private Attributes
DoubleProperty	m_maxDistance
	property, see LArG4GenShowerLib::LArG4GenShowerLib More...
DoubleProperty	m_maxRadius
	property, see LArG4GenShowerLib::LArG4GenShowerLib More...
DoubleProperty	m_minEnergy
	property, see LArG4GenShowerLib::LArG4GenShowerLib More...
DoubleProperty	m_containmentEnergy
	property, see LArG4GenShowerLib::LArG4GenShowerLib More...
DoubleProperty	m_energyFraction
	property, see LArG4GenShowerLib::LArG4GenShowerLib More...
StringArrayProperty	m_lib_struct_files
libMap	m_libraries
	pointer to shower library More...
libMap	m_libraries_by_filename
int	m_stat_numshowers
int	m_stat_valid
int	m_stat_invalid
int	m_stat_nolib
std::map< ShowerLib::IShowerLib *, int >	m_stat_lib_saved
std::map< ShowerLib::IShowerLib *, int >	m_stat_lib_notsaved
std::string	m_physicslist_name
DataObjIDColl	m_extendedExtraObjects
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

Class for shower library generation algorithm.

Create shower library using geant hits

Author

Wolfgang Ehrenfeld, University of Hamburg, Germany

Sasha Glazov, DESY Hamburg, Germany

Version

Id

LArG4GenShowerLib.h 711210 2015-11-27 15:56:00Z jchapman

Definition at line 40 of file LArG4GenShowerLib.h.

Member Typedef Documentation

libMap

typedef std::map<std::string, ShowerLib::IShowerLib*> LArG4GenShowerLib::libMap

private

Definition at line 80 of file LArG4GenShowerLib.h.

StoreGateSvc_t

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

privateinherited

Definition at line 388 of file AthCommonDataStore.h.

Constructor & Destructor Documentation

LArG4GenShowerLib()

LArG4GenShowerLib::LArG4GenShowerLib	(	const std::string &	name,
		ISvcLocator *	pSvcLocator
	)

Definition at line 84 of file LArG4GenShowerLib.cxx.

  : AthAlgorithm(name, pSvcLocator), m_stat_numshowers(0), m_stat_valid(0), m_stat_invalid(0), m_stat_nolib(0)
{
  declareProperty("MaxDistance",   m_maxDistance = 50000.,
                  "max distance squared after which the hits will be truncated");
  declareProperty("MinEnergy",   m_minEnergy = .99,
                  "energy border, that truncation won't cross");
  declareProperty("MaxRadius",        m_maxRadius = 25.,
                  "maximal radius squared until two hits will be combined");
  declareProperty("ContainmentEnergy",        m_containmentEnergy = 0.95,
                  "energy fraction that will be inside containment borders");
  declareProperty("LibStructFiles",   m_lib_struct_files,
                  "List of files to read library structures from");
  declareProperty("EnergyFraction",   m_energyFraction = .02,
                  "the allowed amount of energy that can be deposited outside calorimeter region ");
  declareProperty("PhysicsList", m_physicslist_name = "FTFP_BERT", "Geant4 PhysicsList used in the simulation");
}

Member Function Documentation

addingTagsToLibrary()

void LArG4GenShowerLib::addingTagsToLibrary ( )

private

adding tag information (release, detector description, ...) to library comment

Definition at line 427 of file LArG4GenShowerLib.cxx.

{
  char* atlasProject = getenv("AtlasProject");
  char* atlasVersion = getenv("AtlasVersion");
  std::string atlasReleaseTag = (atlasProject? std::string(atlasProject)+std::string("-") : std::string("Unknown-"));
  if(atlasVersion) {
    atlasReleaseTag += std::string(atlasVersion);
  }
  else {
    atlasReleaseTag += std::string("Unknown");
  }
 
  libMap::iterator itr;
  for (itr = m_libraries.begin();itr != m_libraries.end();++itr){
    // release
    (*itr).second->release(atlasReleaseTag);
 
    // get geometry version
    ServiceHandle<IGeoModelSvc> geoModelSvc("GeoModelSvc", this->name());
    (*itr).second->geometry(geoModelSvc->atlasVersion());
 
    // get Physics list
    (*itr).second->physicsList(m_physicslist_name);
 
    // get geant4 version and strip off CVS Name tag
    std::string g4Version = G4Version;
    size_t pos = g4Version.find("$Name: ");
    if (pos != std::string::npos) {
      g4Version.erase(pos, 7);
    }
    pos = g4Version.find(" $");
    if (pos != std::string::npos) {
      g4Version.erase(pos, 2);
    }
    (*itr).second->geantVersion(g4Version);
  }
}

calculateMoments()

void LArG4GenShowerLib::calculateMoments ( const ShowerLib::StepInfoCollection &  eventSteps, double &  weights, double &  xavfra, double &  yavfra, double &  ravfra  )

private

calculate moments from StepInfoCollection

Definition at line 505 of file LArG4GenShowerLib.cxx.

{
  double escal(0.);
  double xav(0.), yav(0.), xav2(0.), yav2(0.);
 
  for (ShowerLib::StepInfoCollection::const_iterator i=eventSteps.begin();i<eventSteps.end(); ++i) {
    escal += (*i)->dep();
    xav   += (*i)->x()*(*i)->dep();
    yav   += (*i)->y()*(*i)->dep();
    xav2  += (*i)->x()*(*i)->x()*(*i)->dep();
    yav2  += (*i)->y()*(*i)->y()*(*i)->dep();
    //    std::cout << "TRUNC EXYZ"<<(*i)->x()<<" "<<(*i)->y()<<" "<<(*i)->z()
    //      <<" "<<(*i)->dep()<<" "<<(*i)->code<<std::endl;
  }
 
  // Center of gravity:
  const double inv_escal = 1. / escal;
  xavfra = xav*inv_escal;
  yavfra = yav*inv_escal;
  // Second momentum:
  ravfra = std::sqrt(std::abs((xav2*inv_escal-xavfra*xavfra) +
                              (yav2*inv_escal-yavfra*yavfra)));
  // energy is used as weights
  weights = escal;
 
  return;
}

clusterize()

void LArG4GenShowerLib::clusterize ( ShowerLib::StepInfoList *  stepinfo )

private

Definition at line 342 of file LArG4GenShowerLib.cxx.

{
  typedef std::multimap<double,Dist> distMap;
 
  distMap distances;
 
  //fill the map
  for (ShowerLib::StepInfoList::iterator i_h1(stepinfo->begin()); i_h1 != stepinfo->end(); ++i_h1) {
    //iterate only the upper triangle of N*N matrix, since we do not want to create every distance twice
    for (ShowerLib::StepInfoList::reverse_iterator i_h2(stepinfo->rbegin()); (*i_h2) != (*i_h1); ++i_h2) {
      distances.insert(distMap::value_type((*i_h1)->diff2(**i_h2),Dist(*i_h1, *i_h2)));
    }
  }
 
  const double rmin = m_maxRadius;
 
  int cursize = stepinfo->size();
 
 
  while (cursize > 1) {
    distMap::iterator iter = distances.begin();
    while ( !( (*iter).second.h1->valid() && (*iter).second.h2->valid() ) ) {
      distances.erase(iter++); //find the first (i.e. closest) valid pair of hits, removing invalid
    }
 
    if ((*iter).first >  rmin ) break; //if the closest already far enough - stop
 
    ShowerLib::StepInfo* mergedHit = (*iter).second.merge(); //merge two closest hits
 
    for (ShowerLib::StepInfoList::iterator i_h1(stepinfo->begin()); i_h1 != stepinfo->end(); ++i_h1) {
      if ((*i_h1)->valid()){ //only for valid hits
        distances.insert(distMap::value_type((*i_h1)->diff2(*mergedHit),Dist(*i_h1, mergedHit))); //calculate and store distances
      }
    }
 
    stepinfo->push_back(mergedHit); //store the merged hit
    cursize--; //merging == -1 valid hit
  }
 
  // remove invalid
  for (ShowerLib::StepInfoList::iterator i = stepinfo->begin(); i != stepinfo->end();) {
    if ((*i)->valid()) {
      ++i;
    } else {
      delete (*i);
      i = stepinfo->erase(i);
    }
  }
 
}

copyStepInfo()

ShowerLib::StepInfoList * LArG4GenShowerLib::copyStepInfo ( const ShowerLib::StepInfoCollection *  stepinfo )

private

Definition at line 300 of file LArG4GenShowerLib.cxx.

{
  ShowerLib::StepInfoList *rez = new ShowerLib::StepInfoList();
  ShowerLib::StepInfo *copy = nullptr;
 
  for (ShowerLib::StepInfoCollection::const_iterator iter = stepinfo->begin(); iter!=stepinfo->end(); ++iter) {
    copy = new ShowerLib::StepInfo(*(*iter));
    rez->push_back(copy);
  }
 
  return rez;
}

copyStepInfoZeroCleanup()

ShowerLib::StepInfoList * LArG4GenShowerLib::copyStepInfoZeroCleanup ( const ShowerLib::StepInfoCollection *  stepinfo )

private

Definition at line 313 of file LArG4GenShowerLib.cxx.

{
  ShowerLib::StepInfoList *rez = new ShowerLib::StepInfoList();
 
  const double dsame = 1.; // 1mm^2
  ShowerLib::StepInfo *i1 = nullptr;
  ShowerLib::StepInfo *i2 = nullptr;
 
  for (ShowerLib::StepInfoCollection::const_iterator i = stepinfo->begin(); i!=stepinfo->end(); ++i) {
    if(i1 == nullptr) {
      i1 = new ShowerLib::StepInfo(*(*i));
      rez->push_back(i1);
    }
    else {
      i2 = new ShowerLib::StepInfo(*(*i));
      // if distance below cut off, combined and erase
      if ( (i1)->diff2(*i2) < dsame) {
        *i1 += *i2;
        delete i2;
      } else {
        rez->push_back(i2);
        i1 = i2;
      }
    }
  }
 
  return rez;
}

declareGaudiProperty() [1/4]

Gaudi::Details::PropertyBase& AthCommonDataStore< AthCommonMsg< Algorithm > >::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< Algorithm > >::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< Algorithm > >::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< Algorithm > >::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< Algorithm > >::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< Algorithm > >::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< Algorithm > >::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< Algorithm > >::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< Algorithm > >::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< Algorithm > >::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< Algorithm > >::detStore ( ) const

inlineinherited

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

Definition at line 95 of file AthCommonDataStore.h.

{ return m_detStore; }

evtStore() [1/2]

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

inlineinherited

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

Definition at line 85 of file AthCommonDataStore.h.

{ return m_evtStore; }

evtStore() [2/2]

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

inlineinherited

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

Definition at line 90 of file AthCommonDataStore.h.

{ return m_evtStore; }

execute()

StatusCode LArG4GenShowerLib::execute

(

)

Definition at line 144 of file LArG4GenShowerLib.cxx.

{
  const ShowerLib::StepInfoCollection* eventStepsES = getStepInfo();
 
  auto theParticle = getParticleFromMC();
  if (!theParticle) {
    ATH_MSG_ERROR ( "Couldn't get truth particle" );
    return StatusCode::FAILURE;
  }
 
  m_stat_numshowers += 1;
 
  double fraction = eventStepsES->invalid_energy/theParticle->momentum().e();
  if (fraction > m_energyFraction) {
    ATH_MSG_WARNING ( "Shower deposited too much energy outside the calorimeter region (" << (int)(fraction*100) << "%), ignoring" );
    m_stat_invalid += 1;
    return StatusCode::SUCCESS;
  }
  //otherwise shower is valid. even if it ultimately wont go to any library
  m_stat_valid += 1;
 
  ShowerLib::StepInfoList* eventSteps;
 
  if (eventStepsES->size()>500) {
    eventSteps = copyStepInfoZeroCleanup(eventStepsES);
  } else {
    eventSteps = copyStepInfo(eventStepsES);
  }
 
  double etot = 0.;
  for (ShowerLib::StepInfoList::const_iterator iter = eventSteps->begin();iter != eventSteps->end();++iter) {
    etot += (*iter)->dep();
  }
 
  std::stringstream location;
 
  location << eventStepsES->detector << "/" << theParticle->pdg_id();
 
  if (m_libraries.find(location.str()) == m_libraries.end()) {
    ATH_MSG_WARNING ( "No library structure for " << location.str() );
 
    for (ShowerLib::StepInfoList::iterator i(eventSteps->begin());i != eventSteps->end(); ++i) {
      delete (*i);
    }
    delete eventSteps;
    m_stat_nolib += 1;
    return StatusCode::SUCCESS; //not really an error, just lost time
  }
 
  Hep3Vector origin(theParticle->production_vertex()->position().x(),
                    theParticle->production_vertex()->position().y(),
                    theParticle->production_vertex()->position().z());
 
  // Also save direction vector. By default shower lib is created in
  // the direction of the input particle.
  Hep3Vector DirectionShower(theParticle->momentum().px(),
                             theParticle->momentum().py(),
                             theParticle->momentum().pz());
  DirectionShower /= DirectionShower.mag();
  Hep3Vector OrthoShower = DirectionShower.orthogonal();
  OrthoShower /= OrthoShower.mag();
  Hep3Vector CrossShower = DirectionShower.cross(OrthoShower);
  CrossShower /= CrossShower.mag();
 
  // Rotate the hit such that Z direction is along the shower direction
  for (ShowerLib::StepInfoList::iterator i(eventSteps->begin()); i != eventSteps->end(); ++i) {
    Hep3Vector HitDir = (*i)->position() - origin;
 
    (*i)->setX(HitDir*OrthoShower);
    (*i)->setY(HitDir*CrossShower);
    (*i)->setZ(HitDir*DirectionShower);
  }
 
  ATH_MSG_VERBOSE ( "Size of input shower: " << eventSteps->size() );
 
  clusterize(eventSteps);
 
  ATH_MSG_VERBOSE ( "Size after clusterization: " << eventSteps->size() );
 
  if (eventSteps->size() > 10) {
    truncate(eventSteps);
    ATH_MSG_VERBOSE ( "Size after truncation: " << eventSteps->size() );
  }
 
  double maxZ = 0, maxR = 0;
  double edep = 0.;
 
  double containmentEnergy = m_containmentEnergy * etot;
 
  // sort the hits by R, make the border where 95% is deposited
  eventSteps->sort(stepInfoDistCompare(stepInfoDistCompare::R));
  ShowerLib::StepInfoList::const_iterator iter = eventSteps->begin();
  for (;(iter != eventSteps->end()) && (edep < containmentEnergy);++iter) {
    edep += (*iter)->dep();
    maxR = (*iter)->position().r();
  }
 
  edep = 0.;
  // sort the hits by Z, make the border where 95% is deposited
  eventSteps->sort(stepInfoDistCompare(stepInfoDistCompare::Z));
  iter = eventSteps->begin();
  for (;(iter != eventSteps->end()) && (edep < containmentEnergy);++iter) {
    edep += (*iter)->dep();
    maxZ = (*iter)->position().z();
  }
 
  ShowerLib::Shower* shower = new ShowerLib::Shower();
  shower->setZSize(maxZ);
  shower->setRSize(maxR);
 
  for (ShowerLib::StepInfoList::iterator i(eventSteps->begin());i != eventSteps->end(); ++i) {
 
      shower->push_back(new ShowerLib::ShowerEnergySpot(G4ThreeVector((*i)->x(), (*i)->y(), (*i)->z()),(*i)->dep(),(*i)->time()));
 
    delete (*i);
  }
 
  delete eventSteps;
 
  if ((*m_libraries.find(location.str())).second->storeShower(theParticle, shower)) {
    m_stat_lib_saved[(*m_libraries.find(location.str())).second] += 1;
  } else {
    ATH_MSG_WARNING ( "Wasn't able to store shower (" << location.str() << ")" );
    m_stat_lib_notsaved[(*m_libraries.find(location.str())).second] += 1;
  }
 
  ATH_MSG_VERBOSE ( "Successfully finished" );
 
  return StatusCode::SUCCESS;
}

extraDeps_update_handler()

void AthCommonDataStore< AthCommonMsg< Algorithm > >::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

extraOutputDeps()

const DataObjIDColl & AthAlgorithm::extraOutputDeps ( ) const

overridevirtualinherited

Return the list of extra output dependencies.

This list is extended to include symlinks implied by inheritance relations.

Definition at line 50 of file AthAlgorithm.cxx.

{
  // If we didn't find any symlinks to add, just return the collection
  // from the base class.  Otherwise, return the extended collection.
  if (!m_extendedExtraObjects.empty()) {
    return m_extendedExtraObjects;
  }
  return Algorithm::extraOutputDeps();
}

finalize()

StatusCode LArG4GenShowerLib::finalize

(

)

Definition at line 466 of file LArG4GenShowerLib.cxx.

{
 
  // add condition tags
  addingTagsToLibrary();
 
  libMap::iterator itr;
  for (itr = m_libraries_by_filename.begin();itr != m_libraries_by_filename.end();++itr){
    ATH_MSG_DEBUG ( "Writing shower library to file " << (*itr).first );
 
    TFile libr((*itr).first.c_str(),"RECREATE");
 
    if (!(*itr).second->writeToROOT(&libr)) {
      ATH_MSG_ERROR ( "Wasn't able to write " << (*itr).first << ". Probably empty lib." );
    }
 
    libr.Close();
  }
  if (m_stat_numshowers> 0) {
    ATH_MSG_DEBUG ( "********Statistics for GenShowerLib********" );
    ATH_MSG_DEBUG ( "Total number of showers: " << m_stat_numshowers
                   << ", valid: "<< m_stat_valid << " (" << (m_stat_valid*100)/m_stat_numshowers << "%)"
                   << ", invalid: " << m_stat_invalid << " (" << (m_stat_invalid*100)/m_stat_numshowers << "%)" );
    for (itr = m_libraries.begin();itr != m_libraries.end();++itr){
      ATH_MSG_DEBUG ( "*******************************************" );
      std::stringstream ss((*itr).second->statistics());
      for(std::string line; std::getline(ss,line);)
        ATH_MSG_DEBUG ( line );
      ATH_MSG_DEBUG ( "Saved: " << m_stat_lib_saved[(*itr).second] << " Rejected: " << m_stat_lib_notsaved[(*itr).second] );
    }
    ATH_MSG_DEBUG ( "*******************************************" );
    ATH_MSG_DEBUG ( "Showers with no corresponding library: " << m_stat_nolib << " (" << (m_stat_nolib*100)/m_stat_numshowers << "%)" );
  }
  ATH_MSG_DEBUG ( "Finalized." );
 
  return StatusCode::SUCCESS;
}

getParticleFromMC()

HepMC::ConstGenParticlePtr LArG4GenShowerLib::getParticleFromMC ( )

private

return first MC truth particle for event

Definition at line 275 of file LArG4GenShowerLib.cxx.

{
  const McEventCollection* mcEvent;
  if (evtStore()->retrieve(mcEvent,"GEN_EVENT").isFailure()) return nullptr;
 
  // Return the last particle of the event.
  if (mcEvent)
#ifdef HEPMC3
    return !mcEvent->at(0)->particles().empty() ? mcEvent->at(0)->particles().back() : nullptr;
#else
    return ( * (* mcEvent->begin())->particles_end());
#endif
 
  return nullptr;
}

getStepInfo()

const ShowerLib::StepInfoCollection * LArG4GenShowerLib::getStepInfo ( )

private

Definition at line 291 of file LArG4GenShowerLib.cxx.

{
  const ShowerLib::StepInfoCollection* eventStepsES;
  if (evtStore()->retrieve(eventStepsES, "EventSteps").isFailure()) return nullptr;
 
  return eventStepsES;
 
}

initialize()

StatusCode LArG4GenShowerLib::initialize

(

)

Definition at line 103 of file LArG4GenShowerLib.cxx.

{
  ATH_MSG_DEBUG ( "Initializing" );
 
  ShowerLib::IShowerLib* library = nullptr;
 
  std::vector< std::string >::const_iterator nameiter;
 
  ATH_MSG_DEBUG ( "Starting struct files iteration" );
  for (nameiter = m_lib_struct_files.value().begin(); nameiter != m_lib_struct_files.value().end(); ++nameiter ) {
 
    ATH_MSG_DEBUG ( "Struct file: " << (*nameiter) );
 
    library = ShowerLib::iterateStruct(*nameiter);
 
    if (library == nullptr) {
      ATH_MSG_WARNING ( "Library structure file " << (*nameiter) << " doesn't describe a valid library" );
      continue;
    }
 
    std::stringstream location;
    location << library->detector() << "/" << library->particle_id();
    m_libraries[location.str()] = library;
    m_libraries_by_filename[(*nameiter)+".root"] = library;
    m_stat_lib_saved[library]=0;
    m_stat_lib_notsaved[library]=0;
 
    ATH_MSG_VERBOSE ( "Library named " << (*nameiter)+".root" << "is stored at the location " << location.str() );
 
  }
 
  if (m_libraries.empty()) {
    ATH_MSG_ERROR ( "No valid library structure files found. Further execution is pointless." );
    return StatusCode::FAILURE;
  }
 
  ATH_MSG_DEBUG ( "LArG4GenShowerLib " << this->name() << " initialized" );
 
  return StatusCode::SUCCESS;
}

inputHandles()

virtual std::vector<Gaudi::DataHandle*> AthCommonDataStore< AthCommonMsg< Algorithm > >::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< Algorithm >::msg ( ) const

inlineinherited

Definition at line 24 of file AthCommonMsg.h.

                                {
    return this->msgStream();
  }

msg() [2/2]

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

inlineinherited

Definition at line 27 of file AthCommonMsg.h.

                                                  {
    return this->msgStream(lvl);
  }

msgLvl()

bool AthCommonMsg< Algorithm >::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< Algorithm > >::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< Algorithm > >::renounce ( T &  h )

inlineprotectedinherited

Definition at line 380 of file AthCommonDataStore.h.

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

renounceArray()

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

inlineprotectedinherited

remove all handles from I/O resolution

Definition at line 364 of file AthCommonDataStore.h.

                                                          {
    handlesArray.renounce();
  }

sysInitialize()

StatusCode AthAlgorithm::sysInitialize ( )

overridevirtualinherited

Override sysInitialize.

Override sysInitialize from the base class.

Loop through all output handles, and if they're WriteCondHandles, automatically register them and this Algorithm with the CondSvc

Scan through all outputHandles, and if they're WriteCondHandles, register them with the CondSvc

Reimplemented from AthCommonDataStore< AthCommonMsg< Algorithm > >.

Reimplemented in AthAnalysisAlgorithm, AthFilterAlgorithm, PyAthena::Alg, and AthHistogramAlgorithm.

Definition at line 66 of file AthAlgorithm.cxx.

                                       {
  StatusCode sc=AthCommonDataStore<AthCommonMsg<Algorithm>>::sysInitialize();
 
  if (sc.isFailure()) {
    return sc;
  }
  ServiceHandle<ICondSvc> cs("CondSvc",name());
  for (auto h : outputHandles()) {
    if (h->isCondition() && h->mode() == Gaudi::DataHandle::Writer) {
      // do this inside the loop so we don't create the CondSvc until needed
      if ( cs.retrieve().isFailure() ) {
        ATH_MSG_WARNING("no CondSvc found: won't autoreg WriteCondHandles");
        return StatusCode::SUCCESS;
      }
      if (cs->regHandle(this,*h).isFailure()) {
        sc = StatusCode::FAILURE;
        ATH_MSG_ERROR("unable to register WriteCondHandle " << h->fullKey()
                      << " with CondSvc");
      }
    }
  }
  return sc;
}

sysStart()

virtual StatusCode AthCommonDataStore< AthCommonMsg< Algorithm > >::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.

truncate()

void LArG4GenShowerLib::truncate ( ShowerLib::StepInfoList *  stepinfo )

private

Definition at line 393 of file LArG4GenShowerLib.cxx.

{
  double etot = 0;
  stepinfo->sort(stepInfoDistCompare(stepInfoDistCompare::RHO));
 
  for (ShowerLib::StepInfoList::const_iterator i(stepinfo->begin()); i != stepinfo->end(); ++i) {
    etot += (*i)->dep();
  }
 
  double minEnergy = m_minEnergy * etot;
  double rsum = 0.;
  ShowerLib::StepInfoList::iterator iterCut;
  //we're continuing our cycle if  (we're not at the end yet)  && ( (the current hit is still inside maxDistance radius) || (the total energy still less then the threshold) )
  for (iterCut = stepinfo->begin(); (iterCut != stepinfo->end()) && (((*iterCut)->position().mag2() <  m_maxDistance) || (rsum < minEnergy));++iterCut) {
    rsum += (*iterCut)->dep();
  }
 
  if (iterCut == stepinfo->end()) {
    return;
  }
 
  stepinfo->erase(iterCut,stepinfo->end());
 
  if (rsum == 0) { //WTF?
    return;
  }
 
  const double inv_rsum = 1. / rsum;
  for (iterCut = stepinfo->begin(); iterCut != stepinfo->end(); ++iterCut){
    (*iterCut)->setE((*iterCut)->dep() * etot*inv_rsum);
  }
 
}

updateVHKA()

void AthCommonDataStore< AthCommonMsg< Algorithm > >::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);
      }
    }
  }

Member Data Documentation

m_containmentEnergy

DoubleProperty LArG4GenShowerLib::m_containmentEnergy

private

property, see LArG4GenShowerLib::LArG4GenShowerLib

Definition at line 75 of file LArG4GenShowerLib.h.

m_detStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< Algorithm > >::m_detStore

privateinherited

Pointer to StoreGate (detector store by default)

Definition at line 393 of file AthCommonDataStore.h.

m_energyFraction

DoubleProperty LArG4GenShowerLib::m_energyFraction

private

property, see LArG4GenShowerLib::LArG4GenShowerLib

Definition at line 76 of file LArG4GenShowerLib.h.

m_evtStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< Algorithm > >::m_evtStore

privateinherited

Pointer to StoreGate (event store by default)

Definition at line 390 of file AthCommonDataStore.h.

m_extendedExtraObjects

DataObjIDColl AthAlgorithm::m_extendedExtraObjects

privateinherited

Definition at line 79 of file AthAlgorithm.h.

m_lib_struct_files

StringArrayProperty LArG4GenShowerLib::m_lib_struct_files

private

Definition at line 78 of file LArG4GenShowerLib.h.

m_libraries

libMap LArG4GenShowerLib::m_libraries

private

pointer to shower library

Definition at line 81 of file LArG4GenShowerLib.h.

m_libraries_by_filename

libMap LArG4GenShowerLib::m_libraries_by_filename

private

Definition at line 82 of file LArG4GenShowerLib.h.

m_maxDistance

DoubleProperty LArG4GenShowerLib::m_maxDistance

private

property, see LArG4GenShowerLib::LArG4GenShowerLib

Definition at line 72 of file LArG4GenShowerLib.h.

m_maxRadius

DoubleProperty LArG4GenShowerLib::m_maxRadius

private

property, see LArG4GenShowerLib::LArG4GenShowerLib

Definition at line 73 of file LArG4GenShowerLib.h.

m_minEnergy

DoubleProperty LArG4GenShowerLib::m_minEnergy

private

property, see LArG4GenShowerLib::LArG4GenShowerLib

Definition at line 74 of file LArG4GenShowerLib.h.

m_physicslist_name

std::string LArG4GenShowerLib::m_physicslist_name

private

Definition at line 92 of file LArG4GenShowerLib.h.

m_stat_invalid

int LArG4GenShowerLib::m_stat_invalid

private

Definition at line 86 of file LArG4GenShowerLib.h.

m_stat_lib_notsaved

std::map<ShowerLib::IShowerLib*, int> LArG4GenShowerLib::m_stat_lib_notsaved

private

Definition at line 90 of file LArG4GenShowerLib.h.

m_stat_lib_saved

std::map<ShowerLib::IShowerLib*, int> LArG4GenShowerLib::m_stat_lib_saved

private

Definition at line 89 of file LArG4GenShowerLib.h.

m_stat_nolib

int LArG4GenShowerLib::m_stat_nolib

private

Definition at line 87 of file LArG4GenShowerLib.h.

m_stat_numshowers

int LArG4GenShowerLib::m_stat_numshowers

private

Definition at line 84 of file LArG4GenShowerLib.h.

m_stat_valid

int LArG4GenShowerLib::m_stat_valid

private

Definition at line 85 of file LArG4GenShowerLib.h.

m_varHandleArraysDeclared

bool AthCommonDataStore< AthCommonMsg< Algorithm > >::m_varHandleArraysDeclared

privateinherited

Definition at line 399 of file AthCommonDataStore.h.

m_vhka

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

privateinherited

Definition at line 398 of file AthCommonDataStore.h.

---

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

• LArG4GenShowerLib.h
• LArG4GenShowerLib.cxx