702 JParser<> zap(
"Example program to calculate multiples rate.");
707 zap[
'D'] =
make_field(D_m) = JRange_t(0.216, 10);
719 catch(
const exception &error) {
720 FATAL(error.what() << endl);
723 gRandom->SetSeed(seed);
731 DEBUG(module << endl);
735 const double R_m = 17.0 * 2.54 * 0.5e-2;
736 const double A =
PI * R_m * R_m;
738 const double wmin = 280.0;
739 const double wmax = 700.0;
742 const double WAVELENGTH_EXPANSION = (wmax-wmin) / (wmin*wmax) * (300.0*600.0)/(600.0-300.0);
748 case +2: enigma =
new JEnigma<+2>(D_m);
break;
749 case 0: enigma =
new JEnigma< 0>(D_m);
break;
750 case -2: enigma =
new JEnigma<-2>(D_m);
break;
755 const double vmin = 1.0 / wmax;
756 const double vmax = 1.0 / wmin;
760 for (
double w = wmin;
w <= wmax;
w += 1.0) {
766 NOTICE(
"Maximal QE " <<
FIXED(5,3) << QEmax << endl);
767 NOTICE(
"Wavelength expansion " <<
FIXED(5,3) << WAVELENGTH_EXPANSION << endl);
768 NOTICE(
"Number of photons per decay " <<
FIXED(5,2) << ng << endl);
772 JManager_t
H1(
new TH1D(
"M[%]", NULL, 100, D_m.getLowerLimit(), D_m.getUpperLimit()));
776 TH1D pmt(
"pmt", NULL, 1000, -1.0, +1.0);
778 for (Int_t i = 1; i != pmt.GetNbinsX(); ++i) {
780 const double dot = pmt.GetBinCenter(i);
791 y = get_angular_acceptance(dot);
795 pmt.SetBinContent(i, y);
806 for (
counter_type event_count = 0; event_count != numberOfEvents; ++event_count) {
808 if (event_count%10000 == 0) {
809 STATUS(
"event: " << setw(10) << event_count <<
'\r');
DEBUG(endl);
812 const JResult&
result = enigma->next();
814 const double D = result.D;
815 const double V = result.V;
820 double W =
A / (4*
PI*(D-R_m)*(D-R_m));
827 double x = gRandom->Rndm();
830 if ((x -= k40_beta_decay .getBranchingRatio()) <= 0.0)
831 y = k40_beta_decay (gRandom->Rndm());
832 else if ((x -= k40_electron_capture.getBranchingRatio()) <= 0.0)
833 y = k40_electron_capture(gRandom->Rndm());
835 const int N = gRandom->Poisson(y * WAVELENGTH_EXPANSION * QE * W * QEmax * focus);
844 const double ct = gRandom->Uniform(-1.0, +1.0);
845 const double phi = gRandom->Uniform(-
PI, +
PI);
847 const double st = sqrt((1.0 - ct) * (1.0 + ct));
855 for (
int i = 0; i !=
N; ++i) {
859 const double v = gRandom->Uniform(vmin, vmax);
860 const double w = 1.0 /
v;
866 for (
size_t pmt = 0; pmt != module.size(); ++pmt) {
872 const double d = pos.getLength();
877 ERROR(
"Distance " << d <<
" < " << D << endl);
897 p = get_angular_acceptance(dot) *
getQE(w);
901 P += pi[pmt] = U * p *
exp(-d/l_abs);
905 ERROR(
"Probability " << P <<
" > " << W << endl);
908 if (W * QEmax * gRandom->Rndm() <
P) {
911 double y = gRandom->Uniform(P);
915 buffer.push_back(pmt);
919 if (!buffer.empty()) {
921 int M = buffer.size();
925 sort(buffer.begin(), buffer.end());
927 M =
distance(buffer.begin(), unique(buffer.begin(), buffer.end()));
933 for (
int i = 2; i <=
M; ++i) {
934 P2[i].put((
double) (buffer.size() -
M) / (
double)
M, V);
942 for (JManager_t::iterator i =
H1.begin(); i !=
H1.end(); ++i) {
943 i->second->Scale(bequerel / (
double) numberOfEvents);
946 for (
size_t M = 2; M != 7; ++
M) {
947 cout <<
"Rate[" << M <<
"] = "
948 <<
FIXED(7,3) << bequerel *
h1[
M].getTotal() / (double) numberOfEvents
950 <<
FIXED(7,3) << bequerel *
h1[
M].getError() / (double) numberOfEvents
954 for (
size_t M = 2; M != 7; ++
M) {
956 cout <<
"P2[" << M <<
"] = " << P2[
M].getMean() << endl;
Utility class to parse command line options.
const double getPhotocathodeArea()
Photo-cathode area 10 inch PMT.
Data structure for a composite optical module.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
double getDot(const JNeutrinoDirection &first, const JNeutrinoDirection &second)
Dot product.
then for HISTOGRAM in h0 h1
Long64_t counter_type
Type definition for counter.
Auxiliary data structure for floating point format specification.
Abstract interface for the generation of points in 3D space.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` typeset -Z 4 STRING JOpera1D -f hydrophone.root
Description of Monte Carlo event generation applications.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
JDirection3D getDirection(const Vec &dir)
Get direction.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
double getQE(const double R, const double mu)
Get QE for given ratio of hit probabilities and expectation value of the number of photo-electrons...
JPosition3D getPosition(const Vec &pos)
Get position.
static const double PI
Mathematical constants.
JDetectorAddressMap & getDetectorAddressMap()
Get detector address map.
double getAbsorptionLength(const double lambda)
Absoption length.
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
int getCount(const T &hit)
Get hit count.
Data structure for position in three dimensions.
static const JPhotocathodeArea2D getPhotocathodeArea2D
Function object for photo-cathode area 3 inch PMT.
then usage $script[input file[working directory[option]]] nWhere option can be N
do echo Generating $dir eval D
source $JPP_DIR setenv csh $JPP_DIR eval JShellParser o a A
const JModule & getModule(const JDetector &detector, const JModuleLocation &location)
find module with a given string and floor number
double getAngularAcceptance(const double x)
Angular acceptence of PMT.