LCDG4CalSDProj.cc

Go to the documentation of this file.

#include "LCDG4CalSDProj.hh"
#include "LCD_DocumentUtil.hh"
#include "LCD_DocGeoUtil.hh"

// Geant4
#include "G4String.hh"


XERCES_CPP_NAMESPACE_USE;

/*
 * Constructor
 *
 * This constructor creates a cal-type SD that models projective
 * geometries with fixed-solid-angle cells.
 * XML parms are loaded from voldb, which points to the <TYPE>_BARREL
 * volume.  Parameters necessary for processing projective SD in 
 * ProcessHits are stored in phiseg and thetaseg arrays
 */
LCDG4CalSDProj::LCDG4CalSDProj(G4String name, const DOM_Element& voldb)
  : LCDG4CalSD(name)
{
  G4cout <<"  Creating projective LCDG4CalSD..."<< G4endl;
  double thseg, phiseg;
  bool ok = LCD_DocGeoUtil::getSegmentation( voldb, &thseg, &phiseg);
  if (!ok) {
    G4cerr <<"***** LCDG4CalSDProj Error: can not get segmentation!!!"
           << G4endl;
    exit(1);
  }

  _thetaseg = (G4int)thseg;
  _phiseg = (G4int)phiseg;
}

/*
 * Destructor
 */
LCDG4CalSDProj::~LCDG4CalSDProj() {
}

//*** Projective cell indices
unsigned int LCDG4CalSDProj::findCell(const G4ThreeVector& pos) {

  // theta,phi are spherical coords from geant hit
  G4double hittheta = pos.theta();
  G4double hitphi = pos.phi();
  if(hitphi < 0.0) hitphi += 360.0*degree;

  G4double dphi   = 360.0*degree/_phiseg;
  G4double dtheta = 180.0*degree/_thetaseg;
  // convenience aliases
  unsigned int& ilayer = index[1];
  unsigned int& iphi   = index[2];
  unsigned int& itheta = index[3];
  // set geometry-dependent indices: index[2] and index[3]
  iphi   = int(hitphi/dphi);
  itheta = int(hittheta/dtheta);
  unsigned int cellIndex = itheta + _thetaseg* (iphi + _phiseg*ilayer);
  return cellIndex;
}

//*** Projective geometry calculations to find cell center
G4ThreeVector LCDG4CalSDProj::getCellCenter(unsigned int cellIndex) {
  if(cellIndex) ; // make compiler happy
  // convenience aliases
  unsigned int& sysno  = index[0];
  unsigned int& ilayer = index[1];
  unsigned int& iphi   = index[2];
  unsigned int& itheta = index[3];

  double dtheta = 180.*degree / _thetaseg;
  double dphi = 360.*degree / _phiseg;
  double theta = (itheta+0.5)*dtheta/radian;  // +0.5 for center of cell
  double phi = (iphi+0.5)*dphi/radian;        // +0.5 for center of cell

  // rcyl is cyl radius for the center of active material in i-th layer
  // rsph is spherical rad for the center of active
  double rcyl=0, zc=0;
  if(sysno%2==0) {
    // ** BARREL **
    rcyl = layer0CenterDepth + ilayer*layerThickness;
    zc = rcyl/tan(theta);
  }
  else {
    // ** ENDCAP **
    zc = layer0CenterDepth + ilayer*layerThickness;
    // Check for north/south
    if(north) zc = -zc;
    rcyl = zc*tan(theta);
    assert(rcyl >= 0.0);
  }


  double xc = rcyl*cos(phi);
  double yc = rcyl*sin(phi);
  G4ThreeVector cellcenter(xc,yc,zc);
  return cellcenter;
}

// Index decoding: index1 is iphi
unsigned int LCDG4CalSDProj::getIndex1(unsigned int cellIndex) const {
  unsigned int iaux = cellIndex - (cellIndex%_thetaseg);
  assert( iaux%_thetaseg == 0 );
  unsigned int iphi = (iaux/_thetaseg) % _phiseg;
  return iphi;
}

// Index decoding: index2 is itheta
unsigned int LCDG4CalSDProj::getIndex2(unsigned int cellIndex) const {
  return cellIndex % _thetaseg;
}

// Index decoding: index3 is ilayer
unsigned int LCDG4CalSDProj::getIndex3(unsigned int cellIndex) const {
  unsigned int imult = _phiseg*_thetaseg;
  unsigned int iaux = cellIndex - (cellIndex % imult);
  unsigned int ilayer = iaux / imult;
  return ilayer;
}

---