Jpp  pmt_effective_area_update_2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JMultiplicityK40.cc File Reference

Example program to calculate multiples rate. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TRandom3.h"
#include "JTools/JRange.hh"
#include "JTools/JWeight.hh"
#include "JTools/JQuantile.hh"
#include "JPhysics/KM3NeT.hh"
#include "JPhysics/JConstants.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JMath/JMathToolkit.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorAddressMapToolkit.hh"
#include "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Detailed Description

Example program to calculate multiples rate.

Author
mdejong

Definition in file JMultiplicityK40.cc.

Function Documentation

int main ( int  argc,
char *  argv[] 
)

Definition at line 681 of file JMultiplicityK40.cc.

682 {
683  using namespace std;
684  using namespace JPP;
685 
686  typedef long long int counter_type;
687 
688  string outputFile;
689  double bequerel;
690  JRange_t D_m;
691  counter_type numberOfEvents;
692  double QE;
693  UInt_t seed;
694  int generator;
695  double focus;
696  int Aefftype;
697  bool exclusive;
698  int debug;
699 
700  try {
701 
702  JParser<> zap("Example program to calculate multiples rate.");
703 
704  zap['o'] = make_field(outputFile) = "k40.root";
705  zap['b'] = make_field(bequerel) = 13.75e3; // [m^-3 s^-1]
706  zap['n'] = make_field(numberOfEvents) = 1e7;
707  zap['D'] = make_field(D_m) = JRange_t(0.216, 10);
708  zap['Q'] = make_field(QE) = 1.0;
709  zap['S'] = make_field(seed) = 0;
710  zap['G'] = make_field(generator) = 0, +2, -2;
711  zap['F'] = make_field(focus) = 1.0;
712  zap['A'] = make_field(Aefftype) = 1, 2, 3;
713  zap['U'] = make_field(exclusive);
714  zap['a'] = make_field(a) = 0.0;
715  zap['d'] = make_field(debug) = 2;
716 
717  zap(argc, argv);
718  }
719  catch(const exception &error) {
720  FATAL(error.what() << endl);
721  }
722 
723  gRandom->SetSeed(seed);
724 
725  using namespace KM3NET;
726  using namespace JPP;
727 
728  const int id = 1;
729  const JModule module = getModule(getDetectorAddressMap(14).get(id), id);
730 
731  DEBUG(module << endl);
732 
733  bequerel /= focus; // correct decay rate for focussing of photons
734 
735  const double R_m = 17.0 * 2.54 * 0.5e-2; // radius 17'' optical module [m]
736  const double A = PI * R_m * R_m; // cross section optical module [m^2]
737 
738  const double wmin = 280.0; // minimal wavelength [nm]
739  const double wmax = 700.0; // maximal wavelength [nm]
740  double ng = 41.0; // average number of photons per decay in given wavelength range
741 
742  const double WAVELENGTH_EXPANSION = (wmax-wmin) / (wmin*wmax) * (300.0*600.0)/(600.0-300.0);
743 
744  JGenerator* enigma = NULL; // distance generation
745 
746  switch (generator) {
747 
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;
751 
752  default: FATAL("No generator type " << generator << endl);
753  }
754 
755  const double vmin = 1.0 / wmax; // wavelength generation
756  const double vmax = 1.0 / wmin; // as lambda^-2
757 
758  double QEmax = 0.0;
759 
760  for (double w = wmin; w <= wmax; w += 1.0) {
761  if (getQE(w) > QEmax) {
762  QEmax = getQE(w);
763  }
764  }
765 
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);
769 
770  typedef JManager<int, TH1D> JManager_t;
771 
772  JManager_t H1(new TH1D("M[%]", NULL, 100, D_m.getLowerLimit(), D_m.getUpperLimit()));
773 
774  H1->Sumw2();
775 
776  TH1D pmt("pmt", NULL, 1000, -1.0, +1.0);
777 
778  for (Int_t i = 1; i != pmt.GetNbinsX(); ++i) {
779 
780  const double dot = pmt.GetBinCenter(i);
781 
782  double y = 0.0;
783 
784  switch (Aefftype) {
785 
786  case 1:
787  y = getAngularAcceptance(dot);
788  break;
789 
790  case 3:
791  y = get_angular_acceptance(dot);
792  break;
793  }
794 
795  pmt.SetBinContent(i, y);
796  }
797 
798 
801 
802 
803  vector<double> pi(module.size());
804  vector<int> buffer;
805 
806  for (counter_type event_count = 0; event_count != numberOfEvents; ++event_count) {
807 
808  if (event_count%10000 == 0) {
809  STATUS("event: " << setw(10) << event_count << '\r'); DEBUG(endl);
810  }
811 
812  const JResult& result = enigma->next();
813 
814  const double D = result.D;
815  const double V = result.V;
816 
817  // check first photons on cross section of optical module at its face
818  // it is assumed that the light from the K40 decay is purely isotropic
819 
820  double W = A / (4*PI*(D-R_m)*(D-R_m));
821 
822  if (W > 0.5) {
823  W = 0.5;
824  }
825 
826  //
827  double x = gRandom->Rndm(); // decay mode
828  double y = 0.0; // number of photons
829 
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());
834 
835  const int N = gRandom->Poisson(y * WAVELENGTH_EXPANSION * QE * W * QEmax * focus);
836  //
837  /*
838  const int N = gRandom->Poisson(ng * WAVELENGTH_EXPANSION * QE * W * QEmax * focus);
839  */
840  if (N >= 2) {
841 
842  // generate vertex
843 
844  const double ct = gRandom->Uniform(-1.0, +1.0);
845  const double phi = gRandom->Uniform(-PI, +PI);
846 
847  const double st = sqrt((1.0 - ct) * (1.0 + ct));
848 
849  const JPosition3D vertex(D * st * cos(phi),
850  D * st * sin(phi),
851  D * ct);
852 
853  buffer.clear();
854 
855  for (int i = 0; i != N; ++i) {
856 
857  // generate wavelength of photon
858 
859  const double v = gRandom->Uniform(vmin, vmax);
860  const double w = 1.0 / v;
861 
862  const double l_abs = getAbsorptionLength(w);
863 
864  double P = 0.0;
865 
866  for (size_t pmt = 0; pmt != module.size(); ++pmt) {
867 
868  JPosition3D pos(module[pmt].getPosition());
869 
870  pos -= vertex;
871 
872  const double d = pos.getLength();
873 
874  pos /= d; // direction of photon
875 
876  if (d < D - R_m) {
877  ERROR("Distance " << d << " < " << D << endl);
878  }
879 
880  const double dot = getDot(pos, module[pmt].getDirection());
881 
882  const double U = 0.5 * (1.0 - d / sqrt(d*d + getPhotocathodeArea() / PI));
883 
884  double p = 0.0;
885 
886  switch (Aefftype) {
887 
888  case 1:
889  p = getAngularAcceptance(dot) * getQE(w);
890  break;
891 
892  case 2:
894  break;
895 
896  case 3:
897  p = get_angular_acceptance(dot) * getQE(w);
898  break;
899  }
900 
901  P += pi[pmt] = U * p * exp(-d/l_abs);
902  }
903 
904  if (P > W) {
905  ERROR("Probability " << P << " > " << W << endl);
906  }
907 
908  if (W * QEmax * gRandom->Rndm() < P) {
909 
910  int pmt = 0;
911  double y = gRandom->Uniform(P);
912 
913  for (vector<double>::const_iterator i = pi.begin(); i != pi.end() && (y -= *i) > 0.0; ++i, ++pmt) {}
914 
915  buffer.push_back(pmt);
916  }
917  }
918 
919  if (!buffer.empty()) {
920 
921  int M = buffer.size();
922 
923  if (exclusive) {
924 
925  sort(buffer.begin(), buffer.end());
926 
927  M = distance(buffer.begin(), unique(buffer.begin(), buffer.end()));
928  }
929 
930  h1[M].put(V);
931  H1[M]->Fill(D, V);
932 
933  for (int i = 2; i <= M; ++i) {
934  P2[i].put((double) (buffer.size() - M) / (double) M, V);
935  }
936 
937  }
938  }
939  }
940  STATUS(endl);
941 
942  for (JManager_t::iterator i = H1.begin(); i != H1.end(); ++i) {
943  i->second->Scale(bequerel / (double) numberOfEvents);
944  }
945 
946  for (size_t M = 2; M != 7; ++M) {
947  cout << "Rate[" << M << "] = "
948  << FIXED(7,3) << bequerel * h1[M].getTotal() / (double) numberOfEvents
949  << " +/- "
950  << FIXED(7,3) << bequerel * h1[M].getError() / (double) numberOfEvents
951  << " Hz" << endl;
952  }
953 
954  for (size_t M = 2; M != 7; ++M) {
955  if (P2[M].getCount() != 0) {
956  cout << "P2[" << M << "] = " << P2[M].getMean() << endl;
957  }
958  }
959 
960  TFile out(outputFile.c_str(), "recreate");
961 
962  out << H1 << pmt;
963 
964  out.Write();
965  out.Close();
966 }
Utility class to parse command line options.
Definition: JParser.hh:1500
data_type w[N+1][M+1]
Definition: JPolint.hh:741
const double getPhotocathodeArea()
Photo-cathode area 10 inch PMT.
Definition: Antares.hh:41
do $JPP JMEstimator M
Definition: JMEstimator.sh:37
Data structure for a composite optical module.
Definition: JModule.hh:66
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.
Definition: JAstronomy.hh:409
do rm f tmp H1
#define STATUS(A)
Definition: JMessage.hh:63
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
Long64_t counter_type
Type definition for counter.
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
string outputFile
#define KM3NET
Definition: Jpp.hh:32
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.
Definition: JHead.hh:403
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
JDirection3D getDirection(const Vec &dir)
Get direction.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
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...
return result
Definition: JPolint.hh:727
JPosition3D getPosition(const Vec &pos)
Get position.
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
static const double PI
Mathematical constants.
int debug
debug level
Definition: JSirene.cc:63
JDetectorAddressMap & getDetectorAddressMap()
Get detector address map.
#define FATAL(A)
Definition: JMessage.hh:67
double getAbsorptionLength(const double lambda)
Absoption length.
Definition: Antares.hh:53
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
Definition: JMuonPath.sh:47
then JCalibrateToT a
Definition: JTuneHV.sh:116
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.
Definition: JPosition3D.hh:36
data_type v[N+1][M+1]
Definition: JPolint.hh:740
static const JPhotocathodeArea2D getPhotocathodeArea2D
Function object for photo-cathode area 3 inch PMT.
Definition: KM3NeT.hh:5598
then usage $script[input file[working directory[option]]] nWhere option can be N
Definition: JMuonPostfit.sh:36
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
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
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:62
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:84