LCDG4CalSDNonProj.cc

Go to the documentation of this file.

#include "LCDG4CalSDNonProj.hh"
#include "LCDG4CalHit.hh"
#include "LCD_DocumentUtil.hh"
#include "NonProjCalGeometry.hh"
#include "NonProjCalGeomParams.hh"
#include "Util_String.hh"
#include "NPCellID.hh"

// Geant4
#include "G4String.hh"

using std::ofstream;
using std::ios;

XERCES_CPP_NAMESPACE_USE;

/*
 * Constructor
 *
 * This constructor creates a cal-type SD that models non-projective
 * geometries with fixed-size rectangular cells.
 * XML parms are loaded from voldb, which points to the <TYPE>_BARREL
 * volume.  Parameters necessary for processing non-projective SD in 
 * ProcessHits are stored in NonProj geometry classes
 */
LCDG4CalSDNonProj::LCDG4CalSDNonProj(G4String name, const DOM_Element& voldb)
  : LCDG4CalSD(name)
{
  G4cout <<"  Creating non-projective LCDG4CalSD..."<< G4endl;

  // Encapsulate the NonProj parameters
  hcalpar = new NonProjCalGeomParams(voldb);

  // Create the geometry object
  hcal = new NonProjCalGeometry(*hcalpar);

  //------------------------------
  // Write out geometry parameters in output txt file.  Temporary?
  const char* geomfile = getenv("TXTGEOMFILE");

  if (geomfile == "") {
    G4cerr << "geomfilepath is NULL";
  }
  else {
    /* print all cell_info, layering, slice and volume information */
    if(debugLevel) {
      G4cout << "For NonProj geometry parameters, please check:\n";
    }
    G4cout << "NP geom file: " << geomfile << G4endl;
  }

  ofstream* outf_hcal = new ofstream(geomfile,ios::app);
  if( outf_hcal ) {
    *outf_hcal <<"Segmentation:\n\n";
    *outf_hcal<<"CalorimeterCell=hep.lcd.geometry.np.NPHADCalorimeterCell\n\n";
    *outf_hcal << "innerR= "<< hcalpar->GlobalRmin()/cm << G4endl;
    *outf_hcal << "outerR= "<< hcalpar->GlobalRmax()/cm << G4endl;
    *outf_hcal << "outerZ= "<< hcalpar->GlobalZmax()/cm << G4endl;
    *outf_hcal << "cellLenZ= "<< hcalpar->CellDimZ()/cm << G4endl;
    *outf_hcal << "cellLenPhi= "<< hcalpar->CellDimPhi()/cm<<G4endl;
    *outf_hcal << "num_layers= "<< hcalpar->NumberOfLayers() << G4endl;
    *outf_hcal << "preabs_thick= "<< hcalpar->PreAbsThickness()/cm<< G4endl;
    *outf_hcal << "sens_thick= "<< hcalpar->SensorThickness()/cm<< G4endl;
    *outf_hcal << "postabs_thick= "<< hcalpar->PostAbsThickness()/cm<< G4endl;
    *outf_hcal << "\nend;\n";
    outf_hcal->close();
    delete outf_hcal;
  }
}

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

// NonProjective cell indices
unsigned int LCDG4CalSDNonProj::findCell(const G4ThreeVector& pos) {
  NPCellID npid = hcal->getClosestCell(pos);
  unsigned int cellIndex = npid.packIn32Bits();
  return cellIndex;
}

// Nonprojective geometry calculations to find cell center
G4ThreeVector LCDG4CalSDNonProj::getCellCenter(unsigned int cellIndex) {

  // Important: set geometry-dependent indices: index[2] and index[3]
  // convenience aliases
  unsigned int& iphi = index[2];
  unsigned int& iz   = index[3];
  iphi = getIndex1(cellIndex);
  iz   = getIndex2(cellIndex);
  NPCellID npid(cellIndex);
  G4ThreeVector cellcenter(hcal->getCellCenter(npid));
  return cellcenter;
}

//*** Cell indices
// Non Projective cell indices
//    6 bits for layer number
//   12 bits for theta index
//   14 bits for phi index
//     = 32 bits total

// index1 is iphi
unsigned int LCDG4CalSDNonProj::getIndex1(unsigned int cellIndex) const {
  return (cellIndex & 0x3fff);
}

// index2 is iz
unsigned int LCDG4CalSDNonProj::getIndex2(unsigned int cellIndex) const {
//   int izAux = (cellIndex>>14)&0xfff;
//   if(izAux&0x800) izAux |= ~0xfff;        // account for sign bit
//   unsigned int iz = izAux;
  return (cellIndex >> 14) & 0xfff;
}

// index3 is ilayer
unsigned int LCDG4CalSDNonProj::getIndex3(unsigned int cellIndex) const {
  return (cellIndex >> 26) & 0x3f;
}

//   //*** Temporary: check we can recover incoming indices ***
//   int laynoAux = (cellIndex>>26);
//   int izAux = (cellIndex>>14)&0xfff;
//   if(izAux&0x800) izAux |= ~0xfff;        // account for sign bit
//   int iphiAux = cellIndex & 0x3fff;
//   assert( laynoAux == npid.ir );
//   assert( izAux == npid.iz );
//   assert( iphiAux == npid.iphi );
//   //*** end of temporary code

---