/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

DybModifyProperties.cc
Go to the documentation of this file.
00001 
00002 // from STL
00003 #include <algorithm>
00004 #include <cmath>
00005 
00006 // from Gaudi 
00007 #include "GaudiKernel/SvcFactory.h" 
00008 #include "GaudiKernel/MsgStream.h"
00009 #include "GaudiKernel/IDataSelector.h"
00010 #include "GaudiKernel/IConverter.h"
00011 #include "GaudiKernel/IToolSvc.h"
00012 #include "GaudiKernel/IRegistry.h"
00013 #include "GaudiKernel/SmartDataPtr.h"
00014 #include "GaudiKernel/DataObject.h"
00015 #include "GaudiKernel/GenericAddress.h"
00016 #include "GaudiKernel/IDataProviderSvc.h" // djaffe
00017 // from DetDesc 
00018 #include "DetDesc/Material.h" // djaffe
00019 #include "DetDesc/Mixture.h" // djaffe
00020 #include "DetDesc/TabulatedProperty.h" // djaffe
00021 
00022 // from clhep
00023 #include "CLHEP/Units/PhysicalConstants.h" // for hbar*c
00024 
00025 // local
00026 #include "DybModifyProperties.h"
00027 
00028 using namespace CLHEP;
00029 
00030 // ************************************************************
00031 
00032 DybModifyProperties::DybModifyProperties(const std::string& name, ISvcLocator* pSvcLocator):
00033   GaudiAlgorithm(name, pSvcLocator),
00034   m_PropMap ()
00035 {
00036 
00037   // note that PropVector is vector of strings so parameter(s) are in quotes
00038   PropVector v;
00039   v.push_back("ABSLENGTH");
00040   v.push_back("0.007");
00041   m_PropMap[ "/dd/Materials/Water" ] = v;
00042 
00043    // here is an example of how to change the Acrylic index of refraction to a new, constant value of 1.234
00044    // v.push_back("RINDEX");
00045    // v.push_back("NEWCONSTANT");
00046    // v.push_back("1.234")
00047    // m_PropMap[ "/dd/Materials/Acrylic" ] = v;
00048 
00049   // Note that you cannot simultaneously change two properties of one material. 
00050 
00051    declareProperty("ModifyPropertyMap", 
00052                    m_PropMap,
00053                    "Material, Property and parameter(s)");
00054  } 
00055 
00056 // **********************************************************
00057 
00058 StatusCode DybModifyProperties::execute()
00059 {
00060   // only want to do this once
00061   static bool doneit = false;
00062   if ( doneit ) return StatusCode::SUCCESS;
00063   doneit = true;
00064 
00065   debug() << " Hello from DybModifyProperties......" << endreq;
00066   for ( PropMap::iterator it = m_PropMap.begin(); it != m_PropMap.end() ; ++it )
00067     {
00068       const std::string Name = (*it).first;
00069       debug() << " Name of material is " << Name << endreq;
00070 
00071       DataObject* object;
00072       StatusCode sc =  detSvc()->retrieveObject(Name, object);
00073       if ( !object ) {
00074         error() << "Failed to find DataObject with Name " << Name << endreq;
00075         return StatusCode::FAILURE;
00076       }
00077 
00078       Material* mat = dynamic_cast<Material*>( &*object);
00079       if ( !mat ) {
00080         error() << "No material for object made from " << Name << endreq;
00081         return StatusCode::FAILURE;
00082       }
00083 
00084       debug() << "+++++ Original contents of table for " << Name << endreq;
00085       mat->fillStream( debug() ); 
00086 
00087       PropVector pv = (*it).second;
00088       
00089 
00090       for ( SmartRefVector<TabulatedProperty>::const_iterator CIT = mat->tabulatedProperties().begin() ;
00091             CIT != mat->tabulatedProperties().end() ; ++CIT )   
00092         {  if ( pv[0] == (*CIT)->type() )  // change this property
00093             // 
00094             // There are currently two ways to change the property.
00095             // All values in the table can be changed to a new constant 
00096             // or, if the property is ABSLENGTH, the water attenuation algorithm will be applied.
00097             // 
00098             { 
00099               if ( pv[1] == "NEWCONSTANT" ) 
00100                 {
00101 
00102                   double newc = atof(pv[2].c_str()); // should handle units correctly
00103                   
00104                   for ( TabulatedProperty::Table::const_iterator j = ((*CIT)->table()).begin(); j != ((*CIT)->table()).end(); ++j )
00105                     { 
00106                       const double* V = &((*j).second);
00107                       double* val = const_cast<double*> (V);
00108                       *val = newc;
00109                     }
00110                   info() << "For " << Name << " changed all values of " << (*CIT)->type() << " to a constant = " << newc << endreq;
00111 
00112                 }
00113               else
00114                 {
00115                   // Apply the water attenuation algorithm of doc992
00116 
00117                   double a_at_200nm = atof(pv[1].c_str())/cm; // convert string to double
00118                   double defa_at_200nm = 0.007/cm; //  absorbance at 200nm used to construct table in xml file
00119                   
00120                   double lmax = ((*CIT)->table())[0].second; // keep track of max absorption length before, after
00121                   double lmax0 = lmax;
00122                   double wlmax, wlmax0;
00123                   for ( TabulatedProperty::Table::const_iterator j = ((*CIT)->table()).begin();
00124                         j != ((*CIT)->table()).end() ; ++j) // loop over table of (absorption length,photon energy)
00125                     { const double* energy = &((*j).first);
00126                       double wl = twopi * hbarc / (*energy);
00127                       const double* L = &((*j).second);
00128                       double* abslen = const_cast<double*> (L);
00129                       if ( (*abslen)>lmax0 ) { lmax0 = *abslen; wlmax0 = wl;}
00130                       double adef = defa_at_200nm/pow( (wl/(200.*nm)),4 );
00131                       double b = a_at_200nm/pow( (wl/(200.*nm)), 4 );
00132                       
00133                       double anew = 1./(*abslen) - adef + b;
00134                       double ltot = 1./anew;
00135                       if ( ltot > lmax ) { lmax = ltot; wlmax = wl ;}
00136                       *abslen = ltot;
00137                     }
00138                   info() << "For " << Name 
00139                          << " Original max abs. len. "<< lmax0/m << " (m) at " << wlmax0/nm 
00140                          << " (nm) for absorbance(200nm) " << defa_at_200nm*cm 
00141                          << " (1/cm). "<< endreq;
00142                   info() << "For " << Name 
00143                          << " New max abs. len. "<< lmax/m << " (m) at " << wlmax/nm
00144                          << " (nm) for absorbance(200nm) "  << a_at_200nm*cm << " (1/cm)" << endreq;
00145                 }
00146             }
00147         }
00148       debug() << " ++++++ New contents of the table for " << Name << endreq;
00149       mat->fillStream( debug() ); // correct output stream???
00150       
00151     }
00152   return StatusCode::SUCCESS; 
00153 }
00154 
00155 // *****************************************
00156 
00157 StatusCode DybModifyProperties::initialize()
00158 { return StatusCode::SUCCESS; }
00159 
00160 // *****************************************
00161 
00162 StatusCode DybModifyProperties::finalize()
00163 { return StatusCode::SUCCESS; }
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:09:44 for DybAlg by doxygen 1.7.4