0

此问题现已解决 - 修改后的代码如下所示

我在这里有一个问题，我确信只需要少量调整代码，但我似乎无法更正程序。

所以，基本上我想要做的是编写一个 C++ 程序来构造一个 nbin = 20（箱数）的直方图，用于时间间隔 dt（delta t）= 1s 的 10000 个间隔中的盖革计数器的计数数; 假设平均计数率为 5 s^(-1)。为了确定某个时间间隔 deltat 中的计数数，我使用如下形式的 while 语句：

while((t-=tau*log(zscale*double(iran=IM*iran+IC)))<deltat)count++;

作为这个问题的背景知识，我应该提到计数的总数由 n*mu 给出，它与总计数时间 T = n*deltat 成正比。显然，在这个问题中，n 被选择为 10000，deltat 为 1s；给 T = 10000s。

我遇到的问题是我的代码的输出（将在下面显示）只是为元素 0 提供了 10000 个“命中”（对应于时间增量中的 0 个计数），然后当然是 0 个“命中” hist[] 数组的所有其他元素随后。而我期望的输出是泊松分布，其峰值“命中”为 5 个计数（每秒）。

提前感谢您提供的任何帮助，对于我对手头问题的糟糕解释，我深表歉意！我的代码如下所示：

#include <iostream>                                 // Pre-processor directives to include
#include <ctime>                                    //... input/output, time,
#include <fstream>                                  //... file streaming and
#include <cmath>                                    //... mathematical function headers

using namespace std;

int main(void) {

const unsigned IM = 1664525;                    // Integer constants for                
const unsigned IC = 1013904223;                 //... the RNG algorithm
const double zscale = 1.0/0xFFFFFFFF;           // Scaling factor for random double between 0 and 1
const double lambda = 5;                        // Count rate = 5s^-1
const double tau = 1/lambda;                    // Average time tau is inverse of count rate
const int deltat = 1;                           // Time intervals of 1s
const int nbin = 20;                            // Number of bins in histogram
const int nsteps = 1E4;

clock_t start, end;

int count(0);

double t = 0;                               // Time variable declaration

unsigned iran = time(0);                        // Seeds the random-number generator from the system time

int hist[nbin];                                 // Declare array of size nbin for histogram

// Create output stream and open output file

ofstream rout;
rout.open("geigercounterdata.txt");

// Initialise the hist[] array, each element is given the value of zero

for ( int i = 0 ; i < nbin ; i++ )
    hist[i] = 0;

start = clock();

// Construction of histogram using RNG process

for ( int i = 1 ; i <= nsteps ; i++ ) { 

    t = 0;
    count = 0;

    while((t -= tau*log(zscale*double(iran=IM*iran+IC))) < deltat)
        count++;        // Increase count variable by 1

    hist[count]++;          // Increase element "count" of hist array by 1

}

// Print histogram to console window and save to output file

for ( int i = 0 ; i < nbin ; i++ ) {

    cout << i << "\t" << hist[i] << endl;
    rout << i << "\t" << hist[i] << endl;

}

end = clock();

cout << "\nTime taken for process completion = "
     << (end - start)/double(CLOCKS_PER_SEC) 
     << " seconds.\n";

rout.close();

return 1;

}   // End of main() routine
4

1 回答 1

1

我并不完全关注您的 while 循环数学；但是问题确实出在while循环的情况下。我打破了你的while循环如下：

count--;
do 
{
    iran=IM * iran + IC;            //Time generated pseudo-random
    double mulTmp = zscale*iran;    //Pseudo-random double 0 to 1
    double logTmp = log(mulTmp);    //Always negative (see graph of ln(x))
    t -= tau*logTmp;                //Always more than 10^4 as we substract negative
    count++; 
}  while(t < deltat);

从代码中可以明显看出，当您破坏堆时，您总是会遇到count = 0何时t > 1和运行时错误。t < 1

不幸的是，我并没有完全关注你计算背后的数学问题，我不明白为什么会出现泊松分布。对于上述问题，您应该继续解决您的问题（然后为社区分享您的答案），或者为我提供更多数学背景和参考资料，我将使用更正的代码编辑我的答案。如果您决定较早，请记住泊松分布的域是[0, infinity[，因此您需要检查 valecount是否小于 20（或您nbin的）。

于 2013-10-17T21:45:52.950 回答