|
|
|
#include <DreamSD.h>
Public Member Functions | |
| DreamSD (G4String, const DDCompactView &, SensitiveDetectorCatalog &, edm::ParameterSet const &, const SimTrackManager *) | |
| virtual bool | ProcessHits (G4Step *step, G4TouchableHistory *tHistory) |
| virtual uint32_t | setDetUnitId (G4Step *) |
| virtual | ~DreamSD () |
Protected Member Functions | |
| virtual G4bool | getStepInfo (G4Step *aStep) |
| virtual void | initRun () |
Private Types | |
| typedef std::map < G4LogicalVolume *, Doubles > | DimensionMap |
| typedef std::pair< double, double > | Doubles |
Private Member Functions | |
| double | cherenkovDeposit_ (G4Step *aStep) |
| Returns the total energy due to Cherenkov radiation. | |
| const double | crystalLength (G4LogicalVolume *) const |
| const double | crystalWidth (G4LogicalVolume *) const |
| double | curve_LY (G4Step *, int) |
| const double | getAverageNumberOfPhotons_ (const double charge, const double beta, const G4Material *aMaterial, const G4MaterialPropertyVector *rIndex) const |
| Returns average number of photons created by track. | |
| const double | getPhotonEnergyDeposit_ (const G4ParticleMomentum &p, const G4ThreeVector &x, const G4Step *aStep) |
| Returns energy deposit for a given photon. | |
| void | initMap (G4String, const DDCompactView &) |
| bool | setPbWO2MaterialProperties_ (G4Material *aMaterial) |
| Sets material properties at run-time... | |
Private Attributes | |
| double | birk1 |
| double | birk2 |
| double | birk3 |
| std::auto_ptr < G4PhysicsOrderedFreeVector > | chAngleIntegrals_ |
| Table of Cherenkov angle integrals vs photon momentum. | |
| bool | doCherenkov_ |
| G4MaterialPropertiesTable * | materialPropertiesTable |
| int | nphotons_ |
| TTree * | ntuple_ |
| float | px_ [MAXPHOTONS] |
| float | py_ [MAXPHOTONS] |
| float | pz_ [MAXPHOTONS] |
| bool | readBothSide_ |
| int | side |
| double | slopeLY |
| bool | useBirk |
| float | x_ [MAXPHOTONS] |
| DimensionMap | xtalLMap |
| float | y_ [MAXPHOTONS] |
| float | z_ [MAXPHOTONS] |
typedef std::map<G4LogicalVolume*,Doubles> DreamSD::DimensionMap [private] |
typedef std::pair<double,double> DreamSD::Doubles [private] |
| DreamSD::DreamSD | ( | G4String | name, |
| const DDCompactView & | cpv, | ||
| SensitiveDetectorCatalog & | clg, | ||
| edm::ParameterSet const & | p, | ||
| const SimTrackManager * | manager | ||
| ) |
Definition at line 26 of file DreamSD.cc.
References birk1, birk2, birk3, doCherenkov_, g, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), initMap(), edm::Service< T >::isAvailable(), nphotons_, ntuple_, px_, py_, pz_, readBothSide_, slopeLY, useBirk, x_, y_, and z_.
: CaloSD(name, cpv, clg, p, manager) { edm::ParameterSet m_EC = p.getParameter<edm::ParameterSet>("ECalSD"); useBirk= m_EC.getParameter<bool>("UseBirkLaw"); doCherenkov_ = m_EC.getParameter<bool>("doCherenkov"); birk1 = m_EC.getParameter<double>("BirkC1")*(g/(MeV*cm2)); birk2 = m_EC.getParameter<double>("BirkC2"); birk3 = m_EC.getParameter<double>("BirkC3"); slopeLY= m_EC.getParameter<double>("SlopeLightYield"); readBothSide_ = m_EC.getUntrackedParameter<bool>("ReadBothSide", false); edm::LogInfo("EcalSim") << "Constructing a DreamSD with name " << GetName() << "\n" << "DreamSD:: Use of Birks law is set to " << useBirk << " with three constants kB = " << birk1 << ", C1 = " << birk2 << ", C2 = " << birk3 << "\n" << " Slope for Light yield is set to " << slopeLY << "\n" << " Parameterization of Cherenkov is set to " << doCherenkov_ << " and readout both sides is " << readBothSide_; initMap(name,cpv); // Init histogramming edm::Service<TFileService> tfile; if ( !tfile.isAvailable() ) throw cms::Exception("BadConfig") << "TFileService unavailable: " << "please add it to config file"; ntuple_ = tfile->make<TTree>("tree","Cherenkov photons"); if (doCherenkov_) { ntuple_->Branch("nphotons",&nphotons_,"nphotons/I"); ntuple_->Branch("px",px_,"px[nphotons]/F"); ntuple_->Branch("py",py_,"py[nphotons]/F"); ntuple_->Branch("pz",pz_,"pz[nphotons]/F"); ntuple_->Branch("x",x_,"x[nphotons]/F"); ntuple_->Branch("y",y_,"y[nphotons]/F"); ntuple_->Branch("z",z_,"z[nphotons]/F"); } }
| double DreamSD::cherenkovDeposit_ | ( | G4Step * | aStep | ) | [private] |
Returns the total energy due to Cherenkov radiation.
Definition at line 293 of file DreamSD.cc.
References beta, DeDxDiscriminatorTools::charge(), funct::cos(), getAverageNumberOfPhotons_(), getPhotonEnergyDeposit_(), LogDebug, materialPropertiesTable, nphotons_, ntuple_, NULL, phi, px_, py_, pz_, funct::sin(), mathSSE::sqrt(), x_, y_, and z_.
Referenced by getStepInfo().
{
double cherenkovEnergy = 0;
if (!materialPropertiesTable) return cherenkovEnergy;
G4Material* material = aStep->GetTrack()->GetMaterial();
// Retrieve refractive index
const G4MaterialPropertyVector* Rindex = materialPropertiesTable->GetProperty("RINDEX");
if ( Rindex == NULL ) {
edm::LogWarning("EcalSim") << "Couldn't retrieve refractive index";
return cherenkovEnergy;
}
LogDebug("EcalSim") << "Material properties: " << "\n"
<< " Pmin = " << Rindex->GetMinPhotonEnergy()
<< " Pmax = " << Rindex->GetMaxPhotonEnergy();
// Get particle properties
G4StepPoint* pPreStepPoint = aStep->GetPreStepPoint();
G4StepPoint* pPostStepPoint = aStep->GetPostStepPoint();
G4ThreeVector x0 = pPreStepPoint->GetPosition();
G4ThreeVector p0 = aStep->GetDeltaPosition().unit();
const G4DynamicParticle* aParticle = aStep->GetTrack()->GetDynamicParticle();
const double charge = aParticle->GetDefinition()->GetPDGCharge();
// beta is averaged over step
double beta = 0.5*( pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta() );
double BetaInverse = 1.0/beta;
LogDebug("EcalSim") << "Particle properties: " << "\n"
<< " charge = " << charge
<< " beta = " << beta;
// Now get number of photons generated in this step
double meanNumberOfPhotons = getAverageNumberOfPhotons_( charge, beta, material, Rindex );
if ( meanNumberOfPhotons <= 0.0 ) { // Don't do anything
LogDebug("EcalSim") << "Mean number of photons is zero: " << meanNumberOfPhotons
<< ", stopping here";
return cherenkovEnergy;
}
// number of photons is in unit of Geant4...
meanNumberOfPhotons *= aStep->GetStepLength();
// Now get a poisson distribution
int numPhotons = static_cast<int>( G4Poisson(meanNumberOfPhotons) );
//edm::LogVerbatim("EcalSim") << "Number of photons = " << numPhotons;
if ( numPhotons <= 0 ) {
LogDebug("EcalSim") << "Poission number of photons is zero: " << numPhotons
<< ", stopping here";
return cherenkovEnergy;
}
// Material refraction properties
double Pmin = Rindex->GetMinPhotonEnergy();
double Pmax = Rindex->GetMaxPhotonEnergy();
double dp = Pmax - Pmin;
double maxCos = BetaInverse / Rindex->GetMaxProperty();
double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
// Finally: get contribution of each photon
for ( int iPhoton = 0; iPhoton<numPhotons; ++iPhoton ) {
// Determine photon momentum
double randomNumber;
double sampledMomentum, sampledRI;
double cosTheta, sin2Theta;
// sample a momentum (not sure why this is needed!)
do {
randomNumber = G4UniformRand();
sampledMomentum = Pmin + randomNumber * dp;
sampledRI = Rindex->GetProperty(sampledMomentum);
cosTheta = BetaInverse / sampledRI;
sin2Theta = (1.0 - cosTheta)*(1.0 + cosTheta);
randomNumber = G4UniformRand();
} while (randomNumber*maxSin2 > sin2Theta);
// Generate random position of photon on cone surface
// defined by Theta
randomNumber = G4UniformRand();
double phi = twopi*randomNumber;
double sinPhi = sin(phi);
double cosPhi = cos(phi);
// Create photon momentum direction vector
// The momentum direction is still w.r.t. the coordinate system where the primary
// particle direction is aligned with the z axis
double sinTheta = sqrt(sin2Theta);
double px = sinTheta*cosPhi;
double py = sinTheta*sinPhi;
double pz = cosTheta;
G4ThreeVector photonDirection(px, py, pz);
// Rotate momentum direction back to global (crystal) reference system
photonDirection.rotateUz(p0);
// Create photon position and momentum
randomNumber = G4UniformRand();
G4ThreeVector photonPosition = x0 + randomNumber * aStep->GetDeltaPosition();
G4ThreeVector photonMomentum = sampledMomentum*photonDirection;
// Collect energy on APD
cherenkovEnergy += getPhotonEnergyDeposit_( photonMomentum, photonPosition, aStep );
// Ntuple variables
nphotons_ = numPhotons;
px_[iPhoton] = photonMomentum.x();
py_[iPhoton] = photonMomentum.y();
pz_[iPhoton] = photonMomentum.z();
x_[iPhoton] = photonPosition.x();
y_[iPhoton] = photonPosition.y();
z_[iPhoton] = photonPosition.z();
}
// Fill ntuple
ntuple_->Fill();
return cherenkovEnergy;
}
| const double DreamSD::crystalLength | ( | G4LogicalVolume * | lv | ) | const [private] |
Definition at line 270 of file DreamSD.cc.
References xtalLMap.
Referenced by curve_LY(), and getPhotonEnergyDeposit_().
| const double DreamSD::crystalWidth | ( | G4LogicalVolume * | lv | ) | const [private] |
Definition at line 280 of file DreamSD.cc.
References tablePrinter::width, and xtalLMap.
Referenced by getPhotonEnergyDeposit_().
| double DreamSD::curve_LY | ( | G4Step * | aStep, |
| int | flag | ||
| ) | [private] |
Definition at line 239 of file DreamSD.cc.
References crystalLength(), LogDebug, CaloSD::setToLocal(), slopeLY, and CommonMethods::weight().
Referenced by getStepInfo().
{
G4StepPoint* stepPoint = aStep->GetPreStepPoint();
G4LogicalVolume* lv = stepPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
G4String nameVolume= lv->GetName();
double weight = 1.;
G4ThreeVector localPoint = setToLocal(stepPoint->GetPosition(),
stepPoint->GetTouchable());
double crlength = crystalLength(lv);
double localz = localPoint.x();
double dapd = 0.5 * crlength - flag*localz; // Distance from closest APD
if (dapd >= -0.1 || dapd <= crlength+0.1) {
if (dapd <= 100.)
weight = 1.0 + slopeLY - dapd * 0.01 * slopeLY;
} else {
edm::LogWarning("EcalSim") << "DreamSD: light coll curve : wrong distance "
<< "to APD " << dapd << " crlength = "
<< crlength << " crystal name = " << nameVolume
<< " z of localPoint = " << localz
<< " take weight = " << weight;
}
LogDebug("EcalSim") << "DreamSD, light coll curve : " << dapd
<< " crlength = " << crlength
<< " crystal name = " << nameVolume
<< " z of localPoint = " << localz
<< " take weight = " << weight;
return weight;
}
| const double DreamSD::getAverageNumberOfPhotons_ | ( | const double | charge, |
| const double | beta, | ||
| const G4Material * | aMaterial, | ||
| const G4MaterialPropertyVector * | rIndex | ||
| ) | const [private] |
Returns average number of photons created by track.
Definition at line 422 of file DreamSD.cc.
References beta, chAngleIntegrals_, and LogDebug.
Referenced by cherenkovDeposit_().
{
const G4double rFact = 369.81/(eV * cm);
if( beta <= 0.0 ) return 0.0;
double BetaInverse = 1./beta;
// Vectors used in computation of Cerenkov Angle Integral:
// - Refraction Indices for the current material
// - new G4PhysicsOrderedFreeVector allocated to hold CAI's
// Min and Max photon momenta
double Pmin = Rindex->GetMinPhotonEnergy();
double Pmax = Rindex->GetMaxPhotonEnergy();
// Min and Max Refraction Indices
double nMin = Rindex->GetMinProperty();
double nMax = Rindex->GetMaxProperty();
// Max Cerenkov Angle Integral
double CAImax = chAngleIntegrals_->GetMaxValue();
double dp = 0., ge = 0., CAImin = 0.;
// If n(Pmax) < 1/Beta -- no photons generated
if ( nMax < BetaInverse) { }
// otherwise if n(Pmin) >= 1/Beta -- photons generated
else if (nMin > BetaInverse) {
dp = Pmax - Pmin;
ge = CAImax;
}
// If n(Pmin) < 1/Beta, and n(Pmax) >= 1/Beta, then
// we need to find a P such that the value of n(P) == 1/Beta.
// Interpolation is performed by the GetPhotonEnergy() and
// GetProperty() methods of the G4MaterialPropertiesTable and
// the GetValue() method of G4PhysicsVector.
else {
Pmin = Rindex->GetPhotonEnergy(BetaInverse);
dp = Pmax - Pmin;
// need boolean for current implementation of G4PhysicsVector
// ==> being phased out
bool isOutRange;
double CAImin = chAngleIntegrals_->GetValue(Pmin, isOutRange);
ge = CAImax - CAImin;
}
// Calculate number of photons
double numPhotons = rFact * charge/eplus * charge/eplus *
(dp - ge * BetaInverse*BetaInverse);
LogDebug("EcalSim") << "@SUB=getAverageNumberOfPhotons"
<< "CAImin = " << CAImin << "\n"
<< "CAImax = " << CAImax << "\n"
<< "dp = " << dp << ", ge = " << ge << "\n"
<< "numPhotons = " << numPhotons;
return numPhotons;
}
| const double DreamSD::getPhotonEnergyDeposit_ | ( | const G4ParticleMomentum & | p, |
| const G4ThreeVector & | x, | ||
| const G4Step * | aStep | ||
| ) | [private] |
Returns energy deposit for a given photon.
Definition at line 557 of file DreamSD.cc.
References crystalLength(), crystalWidth(), relval_parameters_module::energy, PMTResponse::getEfficiency(), LogDebug, and detailsBasic3DVector::y.
Referenced by cherenkovDeposit_().
{
double energy = 0;
// Crystal dimensions
//edm::LogVerbatim("EcalSim") << p << x;
// 1. Check if this photon goes straight to the APD:
// - assume that APD is at x=xtalLength/2.0
// - extrapolate from x=x0 to x=xtalLength/2.0 using momentum in x-y
G4StepPoint* stepPoint = aStep->GetPreStepPoint();
G4LogicalVolume* lv = stepPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
G4String nameVolume= lv->GetName();
double crlength = crystalLength(lv);
double crwidth = crystalWidth(lv);
double dapd = 0.5 * crlength - x.x(); // Distance from closest APD
double y = p.y()/p.x()*dapd;
LogDebug("EcalSim") << "Distance to APD: " << dapd
<< " - y at APD: " << y;
// Not straight: compute probability
if ( fabs(y)>crwidth*0.5 ) {
}
// 2. Retrieve efficiency for this wavelength (in nm, from MeV)
double waveLength = p.mag()*1.239e8;
energy = p.mag()*PMTResponse::getEfficiency(waveLength);
LogDebug("EcalSim") << "Wavelength: " << waveLength << " - Energy: " << energy;
return energy;
}
| bool DreamSD::getStepInfo | ( | G4Step * | aStep | ) | [protected, virtual] |
Reimplemented from CaloSD.
Definition at line 96 of file DreamSD.cc.
References birk1, birk2, birk3, cherenkovDeposit_(), CaloSD::currentID, curve_LY(), doCherenkov_, CaloSD::edepositEM, CaloSD::edepositHAD, CaloSD::getAttenuation(), TrackInformation::getIDonCaloSurface(), LogDebug, CaloSD::preStepPoint, setDetUnitId(), CaloHitID::setID(), side, CaloSD::theTrack, cond::rpcobgas::time, CaloHitID::unitID(), useBirk, and CommonMethods::weight().
Referenced by ProcessHits().
{
preStepPoint = aStep->GetPreStepPoint();
theTrack = aStep->GetTrack();
G4String nameVolume = preStepPoint->GetPhysicalVolume()->GetName();
// take into account light collection curve for crystals
double weight = 1.;
weight *= curve_LY(aStep, side);
if (useBirk) weight *= getAttenuation(aStep, birk1, birk2, birk3);
edepositEM = aStep->GetTotalEnergyDeposit() * weight;
LogDebug("EcalSim") << "DreamSD:: " << nameVolume << " Side " << side
<<" Light Collection Efficiency " << weight
<< " Weighted Energy Deposit " << edepositEM/MeV
<< " MeV";
// Get cherenkov contribution
if ( doCherenkov_ ) {
edepositHAD = cherenkovDeposit_( aStep );
} else {
edepositHAD = 0;
}
double time = (aStep->GetPostStepPoint()->GetGlobalTime())/nanosecond;
unsigned int unitID= setDetUnitId(aStep);
if (side < 0) unitID++;
TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
int primaryID;
if (trkInfo)
primaryID = trkInfo->getIDonCaloSurface();
else
primaryID = 0;
if (primaryID == 0) {
edm::LogWarning("EcalSim") << "CaloSD: Problem with primaryID **** set by "
<< "force to TkID **** "
<< theTrack->GetTrackID() << " in Volume "
<< preStepPoint->GetTouchable()->GetVolume(0)->GetName();
primaryID = theTrack->GetTrackID();
}
bool flag = (unitID > 0);
G4TouchableHistory* touch =(G4TouchableHistory*)(theTrack->GetTouchable());
if (flag) {
currentID.setID(unitID, time, primaryID, 0);
LogDebug("EcalSim") << "CaloSD:: GetStepInfo for"
<< " PV " << touch->GetVolume(0)->GetName()
<< " PVid = " << touch->GetReplicaNumber(0)
<< " MVid = " << touch->GetReplicaNumber(1)
<< " Unit " << currentID.unitID()
<< " Edeposit = " << edepositEM << " " << edepositHAD;
} else {
LogDebug("EcalSim") << "CaloSD:: GetStepInfo for"
<< " PV " << touch->GetVolume(0)->GetName()
<< " PVid = " << touch->GetReplicaNumber(0)
<< " MVid = " << touch->GetReplicaNumber(1)
<< " Unit " << std::hex << unitID << std::dec
<< " Edeposit = " << edepositEM << " " << edepositHAD;
}
return flag;
}
| void DreamSD::initMap | ( | G4String | sd, |
| const DDCompactView & | cpv | ||
| ) | [private] |
Definition at line 191 of file DreamSD.cc.
References DDFilteredView::addFilter(), DDSpecificsFilter::equals, alcazmumu_cfi::filter, DDFilteredView::firstChild(), i, LogDebug, DDFilteredView::logicalPart(), SensitiveDetector::name, DDName::name(), DDBase< N, C >::name(), DDFilteredView::next(), DDSolid::parameters(), DDSpecificsFilter::setCriteria(), DDSolid::shape(), DDLogicalPart::solid(), python::multivaluedict::sort(), tablePrinter::width, and xtalLMap.
Referenced by DreamSD().
{
G4String attribute = "ReadOutName";
DDSpecificsFilter filter;
DDValue ddv(attribute,sd,0);
filter.setCriteria(ddv,DDSpecificsFilter::equals);
DDFilteredView fv(cpv);
fv.addFilter(filter);
fv.firstChild();
const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
std::vector<G4LogicalVolume *>::const_iterator lvcite;
bool dodet=true;
while (dodet) {
const DDSolid & sol = fv.logicalPart().solid();
std::vector<double> paras(sol.parameters());
G4String name = sol.name().name();
G4LogicalVolume* lv=0;
for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++)
if ((*lvcite)->GetName() == name) {
lv = (*lvcite);
break;
}
LogDebug("EcalSim") << "DreamSD::initMap (for " << sd << "): Solid "
<< name << " Shape " << sol.shape() <<" Parameter 0 = "
<< paras[0] << " Logical Volume " << lv;
double length = 0, width = 0;
// Set length to be the largest size, width the smallest
std::sort( paras.begin(), paras.end() );
length = 2.0*paras.back();
width = 2.0*paras.front();
xtalLMap.insert( std::pair<G4LogicalVolume*,Doubles>(lv,Doubles(length,width)) );
dodet = fv.next();
}
LogDebug("EcalSim") << "DreamSD: Length Table for " << attribute << " = "
<< sd << ":";
DimensionMap::const_iterator ite = xtalLMap.begin();
int i=0;
for (; ite != xtalLMap.end(); ite++, i++) {
G4String name = "Unknown";
if (ite->first != 0) name = (ite->first)->GetName();
LogDebug("EcalSim") << " " << i << " " << ite->first << " " << name
<< " L = " << ite->second.first
<< " W = " << ite->second.second;
}
}
| void DreamSD::initRun | ( | ) | [protected, virtual] |
Reimplemented from CaloSD.
Definition at line 162 of file DreamSD.cc.
References materialPropertiesTable, setPbWO2MaterialProperties_(), and xtalLMap.
{
// Get the material and set properties if needed
DimensionMap::const_iterator ite = xtalLMap.begin();
G4LogicalVolume* lv = (ite->first);
G4Material* material = lv->GetMaterial();
edm::LogInfo("EcalSim") << "DreamSD::initRun: Initializes for material "
<< material->GetName() << " in " << lv->GetName();
materialPropertiesTable = material->GetMaterialPropertiesTable();
if ( !materialPropertiesTable ) {
if ( !setPbWO2MaterialProperties_( material ) ) {
edm::LogWarning("EcalSim") << "Couldn't retrieve material properties table\n"
<< " Material = " << material->GetName();
}
materialPropertiesTable = material->GetMaterialPropertiesTable();
}
}
| bool DreamSD::ProcessHits | ( | G4Step * | step, |
| G4TouchableHistory * | tHistory | ||
| ) | [virtual] |
Reimplemented from CaloSD.
Definition at line 74 of file DreamSD.cc.
References CaloSD::createNewHit(), CaloSD::currentHit, CaloSD::edepositEM, CaloSD::edepositHAD, getStepInfo(), CaloSD::hitExists(), NULL, readBothSide_, and side.
{
if (aStep == NULL) {
return true;
} else {
side = 1;
if (getStepInfo(aStep)) {
if (hitExists() == false && edepositEM+edepositHAD>0.)
currentHit = createNewHit();
if (readBothSide_) {
side = -1;
getStepInfo(aStep);
if (hitExists() == false && edepositEM+edepositHAD>0.)
currentHit = createNewHit();
}
}
}
return true;
}
| uint32_t DreamSD::setDetUnitId | ( | G4Step * | aStep | ) | [virtual] |
Implements CaloSD.
Definition at line 182 of file DreamSD.cc.
References errorMatrix2Lands_multiChannel::id, and LogDebug.
Referenced by getStepInfo().
| bool DreamSD::setPbWO2MaterialProperties_ | ( | G4Material * | aMaterial | ) | [private] |
Sets material properties at run-time...
Definition at line 494 of file DreamSD.cc.
References chAngleIntegrals_, getHLTprescales::index, LogDebug, and asciidump::table.
Referenced by initRun().
{
std::string pbWO2Name("E_PbWO4");
if ( pbWO2Name != aMaterial->GetName() ) { // Wrong material!
edm::LogWarning("EcalSim") << "This is not the right material: "
<< "expecting " << pbWO2Name
<< ", got " << aMaterial->GetName();
return false;
}
G4MaterialPropertiesTable* table = new G4MaterialPropertiesTable();
// Refractive index as a function of photon momentum
// FIXME: Should somehow put that in the configuration
const int nEntries = 14;
double PhotonEnergy[nEntries] = { 1.7712*eV, 1.8368*eV, 1.90745*eV, 1.98375*eV, 2.0664*eV,
2.15625*eV, 2.25426*eV, 2.3616*eV, 2.47968*eV, 2.61019*eV,
2.75521*eV, 2.91728*eV, 3.09961*eV, 3.30625*eV };
double RefractiveIndex[nEntries] = { 2.17728, 2.18025, 2.18357, 2.18753, 2.19285,
2.19813, 2.20441, 2.21337, 2.22328, 2.23619,
2.25203, 2.27381, 2.30282, 2.34666 };
table->AddProperty( "RINDEX", PhotonEnergy, RefractiveIndex, nEntries );
aMaterial->SetMaterialPropertiesTable(table); // FIXME: could this leak? What does G4 do?
// Calculate Cherenkov angle integrals:
// This is an ad-hoc solution (we hold it in the class, not in the material)
chAngleIntegrals_ =
std::auto_ptr<G4PhysicsOrderedFreeVector>( new G4PhysicsOrderedFreeVector() );
int index = 0;
double currentRI = RefractiveIndex[index];
double currentPM = PhotonEnergy[index];
double currentCAI = 0.0;
chAngleIntegrals_->InsertValues(currentPM, currentCAI);
double prevPM = currentPM;
double prevCAI = currentCAI;
double prevRI = currentRI;
while ( ++index < nEntries ) {
currentRI = RefractiveIndex[index];
currentPM = PhotonEnergy[index];
currentCAI = 0.5*(1.0/(prevRI*prevRI) + 1.0/(currentRI*currentRI));
currentCAI = prevCAI + (currentPM - prevPM) * currentCAI;
chAngleIntegrals_->InsertValues(currentPM, currentCAI);
prevPM = currentPM;
prevCAI = currentCAI;
prevRI = currentRI;
}
LogDebug("EcalSim") << "Material properties set for " << aMaterial->GetName();
return true;
}
double DreamSD::birk1 [private] |
Definition at line 58 of file DreamSD.h.
Referenced by DreamSD(), and getStepInfo().
double DreamSD::birk2 [private] |
Definition at line 58 of file DreamSD.h.
Referenced by DreamSD(), and getStepInfo().
double DreamSD::birk3 [private] |
Definition at line 58 of file DreamSD.h.
Referenced by DreamSD(), and getStepInfo().
std::auto_ptr<G4PhysicsOrderedFreeVector> DreamSD::chAngleIntegrals_ [private] |
Table of Cherenkov angle integrals vs photon momentum.
Definition at line 65 of file DreamSD.h.
Referenced by getAverageNumberOfPhotons_(), and setPbWO2MaterialProperties_().
bool DreamSD::doCherenkov_ [private] |
Definition at line 57 of file DreamSD.h.
Referenced by DreamSD(), and getStepInfo().
G4MaterialPropertiesTable* DreamSD::materialPropertiesTable [private] |
Definition at line 66 of file DreamSD.h.
Referenced by cherenkovDeposit_(), and initRun().
int DreamSD::nphotons_ [private] |
Definition at line 69 of file DreamSD.h.
Referenced by cherenkovDeposit_(), and DreamSD().
TTree* DreamSD::ntuple_ [private] |
Definition at line 68 of file DreamSD.h.
Referenced by cherenkovDeposit_(), and DreamSD().
float DreamSD::px_[MAXPHOTONS] [private] |
Definition at line 70 of file DreamSD.h.
Referenced by cherenkovDeposit_(), and DreamSD().
float DreamSD::py_[MAXPHOTONS] [private] |
Definition at line 70 of file DreamSD.h.
Referenced by cherenkovDeposit_(), and DreamSD().
float DreamSD::pz_[MAXPHOTONS] [private] |
Definition at line 70 of file DreamSD.h.
Referenced by cherenkovDeposit_(), and DreamSD().
bool DreamSD::readBothSide_ [private] |
Definition at line 57 of file DreamSD.h.
Referenced by DreamSD(), and ProcessHits().
int DreamSD::side [private] |
Definition at line 62 of file DreamSD.h.
Referenced by getStepInfo(), and ProcessHits().
double DreamSD::slopeLY [private] |
Definition at line 59 of file DreamSD.h.
Referenced by curve_LY(), and DreamSD().
bool DreamSD::useBirk [private] |
Definition at line 57 of file DreamSD.h.
Referenced by DreamSD(), and getStepInfo().
float DreamSD::x_[MAXPHOTONS] [private] |
Definition at line 71 of file DreamSD.h.
Referenced by cherenkovDeposit_(), and DreamSD().
DimensionMap DreamSD::xtalLMap [private] |
Definition at line 60 of file DreamSD.h.
Referenced by crystalLength(), crystalWidth(), initMap(), and initRun().
float DreamSD::y_[MAXPHOTONS] [private] |
Definition at line 71 of file DreamSD.h.
Referenced by cherenkovDeposit_(), and DreamSD().
float DreamSD::z_[MAXPHOTONS] [private] |
Definition at line 71 of file DreamSD.h.
Referenced by cherenkovDeposit_(), and DreamSD().
1.7.3