Dear Hani
Geant 4 can't help you getting the Birks constant. For Geant 4 the
constant is a input value that you have to provide and
which comes from an experimental measurement (literature). Once you have
provide the value you can calculate the reduced using the code snippet
below
that I took from Andrea Dotti. I have also attached a root macro that
shows the suppression curve for a few materials that I found in the
literature.
cheers
Hans
birksc1 = 1.29e2;
birksc2 = 9.59e6;
//
// Birks' Law
// ===========
//
// In the case of Scintillator as active medium, we can
// describe the quenching effects with the Birks' law,
// using the expression and the coefficients taken from
// the paper NIM 80 (1970) 239244 for the organic
// scintillator NE102:
// S*dE/dr
// dL/dr = 
// 1 + C1*(dE/dr)
// with:
// S=1
// C1 = 1.29 x 10^2 g*cm^2*MeV^1
// C2 = 9.59 x 10^6 g^2*cm^4*MeV^2
// These are the same values used by ATLAS TileCal
// and CMS HCAL (and also the default in Geant3).
// You can try different values for these parameters,
// to have an idea on the uncertainties due to them,
// by uncommenting one of the lines below.
// To get the "dE/dr" that appears in the formula,
// which has the dimensions
// [ dE/dr ] = MeV * cm^2 / g
// we have to divide the energy deposit in MeV by the
// product of the step length (in cm) and the density
// of the scintillator:
// rho_scintillator = 1.032 g/cm3 .
// Of course, in our case we use only the denominator
// of the Birks' formula above, because we do not have
// digitization, i.e. we only have energy deposit but
// not conversion to photons.
// Birks should not be applied in the case of gamma
// energy depositions (which happens only for the
// photoelectric process), because in this case the
// step length is related to the photoelectric cross
// section, and not to the distance in which the
// energy is actually deposited, that is what is
// needed in order to determine dE/dx which enters
// in the Birks' law.
// Similarly, for neutron energy depositions (which
// have been introduced in Geant4 8.1 as an effective
// way to describe the elastic nuclei recoils below
// a certain kinetic threshold, set by default to
// 100 keV), the step length is related to the neutron
// elastic crosssection, and not to the real ionization
// range of the recoiled nuclei, which should instead
// be considered for the dE/dx in the Birks' law.
// In the case of neutron energy depositions, the
// most correct way to apply the Birks quench would
// be to eliminate them by setting the kinetic
// threshold to 0.0 (variable "edepLimit" in the
// file "G4HadronElasticPhysics.cc" in
// "geant4/physics_lists/hadronic/Packaging/src/" ),
// so that the recoiled nuclei tracks are always
// generated explicitly. This, of course, costs in
// term of CPU performance.
G4double edep = aStep>GetTotalEnergyDeposit() / CLHEP::MeV;
if (edep == 0.) return false;
G4double rho = aStep>GetPreStepPoint()>GetMaterial()>GetDensity() / (CLHEP::g / CLHEP::cm3);
G4double stepLength = aStep>GetStepLength() / CLHEP::cm;
G4double birksFactor = 1.0;
G4double dedx;
G4double vel;
if (aStep>GetTrack()>GetDefinition() == G4Proton::ProtonDefinition()) {
vel = aStep>GetTrack()>GetVelocity();
}
//Do not apply Birks for gamma deposits!
if (stepLength > 1.0e8)//Check, cut if necessary.
{
dedx = edep / (rho * stepLength); //[MeV*cm^2/g]
birksFactor = 1.0 / (1.0 + birksc1 * dedx + birksc2 * dedx * dedx);
if (aStep>GetTrack()>GetDefinition() == G4Gamma::GammaDefinition()) {
birksFactor = 1.0;
}
if (aStep>GetTrack()>GetDefinition() == G4Neutron::NeutronDefinition()) {
birksFactor = 1.0;
}
}
G4double eobsbirks = edep*birksFactor;
On 12/25/2013 02:43 AM, Hani Negm wrote:
> *** Discussion title: Processes Involving Optical Photons
>
> Dear All;
>
> How can I calculate the Birks attenuation constant for scintillator like
> LaBr3(Ce)?
>
> Hani
>
> 
> Visit this GEANT4 at hypernews.slac.stanford.edu message (to reply or unsubscribe) at:
> http://hypernews.slac.stanford.edu/HyperNews/geant4/get/opticalphotons/311/2.html
#include"TF1.h"
#include"TLegend.h"
#include"TAxis.h"
// dedx in [MeV*cm^2/g]
// values for the scintillator Birks' coefficient C1, with C2 = 0.
//C1_scintillator = 1.31e2; C2_scintillator = 0.0; // Standard Birks law.
//C1_scintillator = 8.35e3; C2_scintillator = 0.0; // Zeus SCSN38, lower limit (35%) IEEE TNS Vol. 39 NO.4 (1992), 511
//C1_scintillator = 1.59e2; C2_scintillator = 0.0; // Pilot B (same paper), upper limit (+21%)
// BGO: kB = 6.5 mum/MeV // NIM A439 (2000) 158166
// with rho =7.13 g/cm^3> C1=4.6345e2 C2=0.0
// CsI(Tl): kB = 1.52 mum/MeV // NIM A439 (2000) 158166
// with rho = 4.51 g/cm^3> C1=6.8552e4 C2=0.0
// GSO(Ce): kB = 5.25 mum/MeV // NIM A439 (2000) 158166
// with rho = 6.7 g/cm^3> C1=3.5175e3 C2=0.0
double birksf(double * x,double *p)
{
double dedx = x[0];
return 1.0 / (1.0 + p[0] * dedx + p[1] * dedx * dedx);
}
void birks()
{
TF1 *bf4 = new TF1("birks4",birksf,0,100,2);
bf4>SetParameters(4.6345e2,0.0);
bf4>SetLineWidth(2);
bf4>SetLineColor(7);
bf4>GetXaxis()>SetTitle("dE/dx [MeV/g cm^{2}]");
bf4>GetYaxis()>SetTitle("Birks Factor");
bf4>Draw();
TF1 *bf = new TF1("birks",birksf,0,100,2);
bf>SetParameters(1.29e2,9.59e6);
bf>SetLineWidth(2);
bf>SetLineColor(3);
bf>Draw("SAME");
TF1 *bf1 = new TF1("birks1",birksf,0,100,2);
bf1>SetParameters(8.35e3,0.0);
bf1>SetLineWidth(2);
bf1>SetLineColor(4);
bf1>Draw("SAME");
TF1 *bf2 = new TF1("birks2",birksf,0,100,2);
bf2>SetParameters(1.31e2,0.0);
bf2>SetLineWidth(2);
bf2>SetLineColor(2);
bf2>Draw("SAME");
TF1 *bf3 = new TF1("birks3",birksf,0,100,2);
bf3>SetParameters(1.59e2,0.0);
bf3>SetLineWidth(2);
bf3>SetLineColor(6);
bf3>Draw("SAME");
TF1 *bf5 = new TF1("birks5",birksf,0,100,2);
bf5>SetParameters(6.8552e4,0.0);
bf5>SetLineWidth(2);
bf5>SetLineColor(21);
bf5>Draw("SAME");
TF1 *bf6 = new TF1("birks6",birksf,0,100,2);
bf6>SetParameters(3.5175e3,0.0);
bf6>SetLineWidth(2);
bf6>SetLineColor(8);
bf6>Draw("SAME");
TLegend *legend = new TLegend(.44, .64, .99, .99);
legend>AddEntry(bf, "CMS/ATLAS");
legend>AddEntry(bf2, "Standard Birks law.");
legend>AddEntry(bf1, "Zeus SCSN38, lower limit (35%) IEEE TNS Vol. 39 NO.4 (1992), 511");
legend>AddEntry(bf3, "Zeus SCSN38, upper limit (+21%) IEEE TNS Vol. 39 NO.4 (1992), 511");
legend>AddEntry(bf4, "BGO, NIM A439 (2000) 158166");
legend>AddEntry(bf5, "CsI(Tl), NIM A439 (2000) 158166");
legend>AddEntry(bf6, "GSO(Ge), NIM A439 (2000) 158166");
legend>Draw();
}
