DybModifyProperties.cc

Go to the documentation of this file.

// from STL
#include <algorithm>
#include <cmath>

// from Gaudi 
#include "GaudiKernel/SvcFactory.h" 
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/IDataSelector.h"
#include "GaudiKernel/IConverter.h"
#include "GaudiKernel/IToolSvc.h"
#include "GaudiKernel/IRegistry.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/DataObject.h"
#include "GaudiKernel/GenericAddress.h"
#include "GaudiKernel/IDataProviderSvc.h" // djaffe
// from DetDesc 
#include "DetDesc/Material.h" // djaffe
#include "DetDesc/Mixture.h" // djaffe
#include "DetDesc/TabulatedProperty.h" // djaffe

// from clhep
#include "CLHEP/Units/PhysicalConstants.h" // for hbar*c

// local
#include "DybModifyProperties.h"

using namespace CLHEP;

// ************************************************************

DybModifyProperties::DybModifyProperties(const std::string& name, ISvcLocator* pSvcLocator):
  GaudiAlgorithm(name, pSvcLocator),
  m_PropMap ()
{

  // note that PropVector is vector of strings so parameter(s) are in quotes
  PropVector v;
  v.push_back("ABSLENGTH");
  v.push_back("0.007");
  m_PropMap[ "/dd/Materials/Water" ] = v;

   // here is an example of how to change the Acrylic index of refraction to a new, constant value of 1.234
   // v.push_back("RINDEX");
   // v.push_back("NEWCONSTANT");
   // v.push_back("1.234")
   // m_PropMap[ "/dd/Materials/Acrylic" ] = v;

  // Note that you cannot simultaneously change two properties of one material. 

   declareProperty("ModifyPropertyMap", 
                   m_PropMap,
                   "Material, Property and parameter(s)");
 } 

// **********************************************************

StatusCode DybModifyProperties::execute()
{
  // only want to do this once
  static bool doneit = false;
  if ( doneit ) return StatusCode::SUCCESS;
  doneit = true;

  debug() << " Hello from DybModifyProperties......" << endreq;
  for ( PropMap::iterator it = m_PropMap.begin(); it != m_PropMap.end() ; ++it )
    {
      const std::string Name = (*it).first;
      debug() << " Name of material is " << Name << endreq;

      DataObject* object;
      StatusCode sc =  detSvc()->retrieveObject(Name, object);
      if ( !object ) {
        error() << "Failed to find DataObject with Name " << Name << endreq;
        return StatusCode::FAILURE;
      }

      Material* mat = dynamic_cast<Material*>( &*object);
      if ( !mat ) {
        error() << "No material for object made from " << Name << endreq;
        return StatusCode::FAILURE;
      }

      debug() << "+++++ Original contents of table for " << Name << endreq;
      mat->fillStream( debug() ); 

      PropVector pv = (*it).second;
      

      for ( SmartRefVector<TabulatedProperty>::const_iterator CIT = mat->tabulatedProperties().begin() ;
            CIT != mat->tabulatedProperties().end() ; ++CIT )   
        {  if ( pv[0] == (*CIT)->type() )  // change this property
            // 
            // There are currently two ways to change the property.
            // All values in the table can be changed to a new constant 
            // or, if the property is ABSLENGTH, the water attenuation algorithm will be applied.
            // 
            { 
              if ( pv[1] == "NEWCONSTANT" ) 
                {

                  double newc = atof(pv[2].c_str()); // should handle units correctly
                  
                  for ( TabulatedProperty::Table::const_iterator j = ((*CIT)->table()).begin(); j != ((*CIT)->table()).end(); ++j )
                    { 
                      const double* V = &((*j).second);
                      double* val = const_cast<double*> (V);
                      *val = newc;
                    }
                  info() << "For " << Name << " changed all values of " << (*CIT)->type() << " to a constant = " << newc << endreq;

                }
              else
                {
                  // Apply the water attenuation algorithm of doc992

                  double a_at_200nm = atof(pv[1].c_str())/cm; // convert string to double
                  double defa_at_200nm = 0.007/cm; //  absorbance at 200nm used to construct table in xml file
                  
                  double lmax = ((*CIT)->table())[0].second; // keep track of max absorption length before, after
                  double lmax0 = lmax;
                  double wlmax, wlmax0;
                  for ( TabulatedProperty::Table::const_iterator j = ((*CIT)->table()).begin();
                        j != ((*CIT)->table()).end() ; ++j) // loop over table of (absorption length,photon energy)
                    { const double* energy = &((*j).first);
                      double wl = twopi * hbarc / (*energy);
                      const double* L = &((*j).second);
                      double* abslen = const_cast<double*> (L);
                      if ( (*abslen)>lmax0 ) { lmax0 = *abslen; wlmax0 = wl;}
                      double adef = defa_at_200nm/pow( (wl/(200.*nm)),4 );
                      double b = a_at_200nm/pow( (wl/(200.*nm)), 4 );
                      
                      double anew = 1./(*abslen) - adef + b;
                      double ltot = 1./anew;
                      if ( ltot > lmax ) { lmax = ltot; wlmax = wl ;}
                      *abslen = ltot;
                    }
                  info() << "For " << Name 
                         << " Original max abs. len. "<< lmax0/m << " (m) at " << wlmax0/nm 
                         << " (nm) for absorbance(200nm) " << defa_at_200nm*cm 
                         << " (1/cm). "<< endreq;
                  info() << "For " << Name 
                         << " New max abs. len. "<< lmax/m << " (m) at " << wlmax/nm
                         << " (nm) for absorbance(200nm) "  << a_at_200nm*cm << " (1/cm)" << endreq;
                }
            }
        }
      debug() << " ++++++ New contents of the table for " << Name << endreq;
      mat->fillStream( debug() ); // correct output stream???
      
    }
  return StatusCode::SUCCESS; 
}

// *****************************************

StatusCode DybModifyProperties::initialize()
{ return StatusCode::SUCCESS; }

// *****************************************

StatusCode DybModifyProperties::finalize()
{ return StatusCode::SUCCESS; }