monte carlo.txt

MonteCarlo仿真的C++源码

Monte Carlo Simulation C++ Source Code
This source code was compiled and run on a Sun Sparc Ultra170 running Solaris
2.5. The GNU C compiler gcc version 2.7 was used by calling g++ so that the C++
code would be recognized. The code (named spectrum.cc) was compiled with the
command:

g++ -o spectrum spectrum.cc

The math.h library must be linked in order for the program to successfully run.
This may require the -lm option on some systems.

The minimum information needed to run the program is input as:

spectrum voltage attenuation photoelectrons

The total supply voltage is given by voltage. The attenuation of the
photomultiplier signal is given by attenuation, and is equal to the factor by
which the signal is to be divided. The average number of photoelectrons ejected
by the photocathode is
given by photoelectrons. More options are available by entering a non-zero
number after the other command line parameters. The user is then prompted for a
threshold setting (default is zero), which is multiplied by ten to get the ADC
channel for the
cutoff (the multiplication by ten is convenient because the data output by the
program was rebinned by a factor of ten before being analyzed). The user is then
prompted for a gain factor (default is zero) to fit the program to individual
8854 tubes
because the have different gains. The gain factor is given in percent divided by
ten by which the gain is to be scaled, for example an input of 10 increases the
simulated gain by 100% and an input of -1 decreases the gain by 10%.

A directory called ``newspec'' must be created as a subdirectory of the
directory in which the program is run before running the program. All output
files are put in this directory. Each run of the program results in 31 files.
The format of the file
names is:

voltage_attenuation_photoelectrons_number

The last part of the file name indicates the number of initial photoelectrons
represented in that particular histogram file. The file ending in 0 is the sum
of all other histograms, the spectrum actually seen if the experiment was run.
The histograms
for specific numbers of initial photoelectrons were created to analyze the
effect of the algorithm used to simulate the low charge events on larger numbers
of photoelectrons and to look at the relations between mean and standard
deviations of the
individual peaks without interference from other peaks. The output of the
program can be analyzed in any program that can import text and make histograms.
The file format is one channel per line (200 channels at 2.5 pC per channel).


#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <sys/types.h>
#include <time.h>
#include <math.h>
#include <String.h>

//DEFINES

//number of bins (in x direction)
#define BINS 200

//maximum number of photoelectrons considered
#define MAX_PE 30

//how many events to simulate
#define NUM_EVENTS 40000

//FUNCTION PROTOTYPES

//return a random integer between lower and upper
int myRandInt(int lower, int upper);

//return a random double between lower and upper
double myRandDouble(double lower, double upper);

//return the value of the poisson distribution with
//average value at x
double myPoisson(int x, double avg);

//return value of gauss distribution with mean,
//standard deviation  at x
double myGauss(double x, double mean, double stdev);

//return x! (return is a double to be able to handle
//large numbers)
double factorial(int x);

//return a random value conforming to a poisson
//distribution
int poissonBin(double avg, int xMax, double yMax);

//return a random value conforming to a gauss
//distribution
double gaussBin(double mean, double stdev, double lc);

//calculates unattenuated mean from voltage
double myMean(double voltage);

//calculates unattenuated stdev from voltage
double myStdev(double voltage);

//calculates low charge events for a gauss with mean
double myLowC(double mean);

//returns the value attenuated by att
double myAtt(double value, double att);

//Converts an integer to a String
String myIntToString(int x);

//Converts a double to string, accurate to tenths place
String myTenthsToString(double x);

//Converts a String to a double
double myStringToDouble(String num);

//for optimization calculates the highest x value to
//guess in the poisson
int maxPoissonX(double avg);

//for optimization calculates the highest y value to
//guess in the poisson
double maxPoissonY(double avg);

//rounding integer routine
int myRound(double x);

//MAIN

/*This is the main part of the program.  There are 3
 required command line parameters and one optional.
 They are [voltage] [attenuation] [average number of
 photoelectrons] [editor mode (non zero value)](opt.)
 If the last parameter is omitted a default value of
 zero is used.
  */

void main(int argc, char *argv[]) {
  if(argc < 4) {
    cerr << "Error: too few parameters\n";
    exit(-1);
  }

  srand48(time(NULL));

  /*make 2 D array to hold all of the info to be written
    to files (spectra) there are MAX_PE+1 different spectra
    (one for the combined spectra) and each has BINS number
    of bins
    */

  int spectra[MAX_PE+1][BINS];
  for(int i=0; i<MAX_PE+1; i++)
    for(int j=0; j<BINS; j++)
      spectra[i][j]=0;

  /*
    gf is gain factor, gain factor scales the average and
    stdev. linearly in an attempt to account for variations
    in the gain of the tubes.
  */

  double voltage, att, avgpe, opt;

  opt=0;

  voltage = myStringToDouble(argv[1]);
  att = myStringToDouble(argv[2]);
  while(att <= 0) {
    cout << "Please enter a positive attenuation factor: ";
    cin >> att;
  }
  avgpe = myStringToDouble(argv[3]);

  if(argc > 4)
    opt=myStringToDouble(argv[4]);

  double thresh=0,gf=0;

  if(opt !=0){
    cout << "\n***Entering Option Editing Mode***\n\n"
         << "Enter the threshold (ADC Channel/10)(default=0):";
    cin >> thresh;
    cout << "Enter the gain factor(default=0):";
    cin >> gf;
  }

  double stdev = myStdev(voltage);
  double mean = myMean(voltage);
  int poissonX = maxPoissonX(avgpe);
  double poissonY = maxPoissonY(avgpe);

  /*adjust mean and stdev*/

  mean = mean + gf*mean;
  stdev = stdev + 0.34*gf*stdev;

  /*get Low Charge Fraction before mean is attenuated*/

  double lc=myLowC(mean);

  /*attenuate mean and stdev*/

  mean=myAtt(mean,att);
  stdev=myAtt(stdev,att);

  for(int i=0; i<NUM_EVENTS;) {

    int numPE = poissonBin(avgpe, poissonX, poissonY);

    double tmpGauss=0.0;
    for(int j=0; j<numPE; j++) {
      tmpGauss += gaussBin(mean,stdev,lc);

    }

    if(tmpGauss>=thresh) {
      int gauss = myRound(tmpGauss);
      if(gauss > BINS-1)
        gauss=BINS-1;
      spectra[numPE][gauss]++;
      i++;
    }
  }


  for(int i=0; i<BINS; i++)
    for(int j=1; j<MAX_PE+1; j++)
      spectra[0][i] += spectra[j][i];


  //output is "newspec/[voltage]_[attenuation]_[avgPE]"
  //all to tenths place

  for(int i=0; i<MAX_PE+1; i++) {
    ofstream file_out;
    String name;
    name = "newspec/";
    name += myTenthsToString(voltage);
    name += "_";
    name += myTenthsToString(att);
    name += "_";
    name += myTenthsToString(avgpe);
    name +="_";
    name += myIntToString(i);
    file_out.open(name);

    for(int j=0; j<BINS; j++)
      file_out << spectra[i][j] << endl;
    file_out.close();
  }

}


int myRandInt(int lower, int upper) {
  int done=0;
  int result;
  while(!done) {
    result = int((upper+1-lower)*drand48() + lower);
    if(result != upper+1)
      done=1;
  }
  return result;
}


double myRandDouble(double lower, double upper) {
  return (upper-lower)*drand48() + lower;
}

double myPoisson(int x, double avg) {
  return pow(avg,x)*exp(-avg)/factorial(x);
}

double myGauss(double x, double mean, double stdev) {
  return exp(-0.5*pow((mean-x)/(stdev),2));
}

double factorial(int x) {
  double result=1;
  for(int i=x; i>1; i--)
    result = result*i;
  return result;
}

int poissonBin(double avg, int xMax, double yMax) {
  int done=0;
  int x;
  while(!done) {
    double y = myRandDouble(0,yMax);
    x = myRandInt(0,xMax);
    if(y < myPoisson(x,avg))
      done=1;
  }
  return x;
}


double gaussBin(double mean, double stdev, double lc) {
  int done=0;
  double x;
  while(!done) {
    double y = myRandDouble(0,1);
    x = myRandDouble(0,mean+6*stdev);

    if((x<mean) && (y<lc))
      done=1;

    if(y < myGauss(x,mean,stdev))
      done=1;
  }
  return x;
}

//fit from data
double myMean(double voltage) {
  return exp(0.00379*voltage-6.45);
}

//fit from data
double myStdev(double voltage) {
  return exp(0.0034*voltage-7.05);
}

//fit to data
double myLowC(double mean) {
  return (0.35 - 0.024*exp(0.0091*mean));
}

double myAtt(double value, double att) {
  return value/att;
}


String myIntToString(int x) {
  String result;

  if(x/10==0) {
    result = char('0'+x);
  }
  else {
    result = myIntToString(x/10);
    result += char('0' + x%10);
  }
  return result;
}


String myTenthsToString(double x) {
  int tmp = int(x*10);
  String result = myIntToString(tmp);
  char last = result[result.length()-1];
  result[result.length()-1]='.';
  result += last;
  return result;
}

double myStringToDouble(String num) {
  int decimalq=0;
  int power=-1;
  double result=0;
  for(int i=0; i<num.length(); i++) {
    if(num[i]=='.')
      decimalq=1;
    else if(!decimalq)
      result=10*result + double(num[i]-'0');
    else {
      result=result + double(num[i]-'0')*pow(10,power);
      power--;
    }
  }
  return result;
}

int maxPoissonX(double avg) {
  double maxY = maxPoissonY(avg);
  int i=0;
  for(i = int(avg)+1; (myPoisson(i,avg) > 0.001*maxY)
   && (i<MAX_PE); i++);
  return i;
}


double maxPoissonY(double avg) {

  int tmp = int(avg);
  double max = 0.0;
  int i=tmp-1;

  if(i < 0)
    i=0;
  for(i; i<=tmp+1; i++)
    if(myPoisson(i,avg) > max)
      max = myPoisson(i,avg);
  return max;
}

int myRound(double x) {
  int result = int(x);
  if( (x - result) > 0.5 )
    result++;
  return result;
}