1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
|
#include <iostream>
#include <fstream> // Relevant libraries and particle header called
#include <cstdlib>
#include <cmath>
#include "detector.h"
#include "/Users/Vito/Documents/Section 4/Section 4/particleclass.h"
using namespace std;
int main()
{
ofstream positronfile("/Users/Vito/Documents/Section 5/Section 5/positron_Ek.txt"),
muonfile("/Users/Vito/Documents/Section 5/Section 5/muon_Ek.txt"),
michelfile("/Users/Vito/Documents/Section 5/Section 5/michel_Ek.txt");
// Initialising particles
particle muon (105.7, 0, 0, 29.78), positron (0.511, 0, 0, 69.80);
// Distribution for perfect detector, 1000 muons and positrons
particle marr[1000], parr[1000];
for (int j = 0; j <= 999; j++) // For-loop that fills all 1000 particle array slots generated
{
threevector m_rand (muon.mag(), random_theta(), random_phi(), 'r'), // Particle position
p_rand (positron.mag(), random_theta(), random_phi(), 'r');
marr[j] = particle (105.7, m_rand); // Assigning the array of particles
parr[j] = particle (0.511, p_rand);
muonfile << marr[j].Ek() <<'\t'<< m_rand.gettheta() <<'\t'<< m_rand.getphi() << endl;
positronfile <<parr[j].Ek()<<'\t'<< p_rand.gettheta()<<'\t'<< p_rand.getphi() << endl;
}
// Michel distribution for 1000 positrons (muon decay)
particle michel_parr[1000];
for (int j = 0; j <= 999; j++)
{
threevector p_michel (michel(), random_theta(), random_phi(), 'r'); // Position, again
michel_parr[j] = particle(0.511, p_michel); // Assiging the particles to the array
michelfile << michel_parr[j].Ek() <<'\t' // Find theta at each array location
<< michel_parr[j].gettheta() <<'\t' // NOTE: still works, *exactly*, if you use
<< michel_parr[j].getphi() << endl; // p_michel rather than michel_parr
}
// Using acceptance
particle accept_marr[1000], accept_parr[1000];
detector measure (40, 140, 'd');
double pcounter = 0.0, mcounter = 0.0; // Counter for pion decay and muon decay
for (int j = 0; j <= 999; j++)
{
cout << measure.getEk(parr[j]) << endl;
if(measure.getEk(parr[j])) // If Ek is in direction limits, returns Ek
{
mcounter += mcounter;
}
//cout << measure.getEk(michel_parr[j]) << endl;
if(measure.getEk(michel_parr[j]) > 0.0) // Muon decay particle energy
{
pcounter += pcounter;
}
}
return 0;
}
// Custom function declarations
double michel() // Creating the michel function
{
double x, y;
do
{
x = 53.0*rand()/RAND_MAX;
y = x/53.0;
}
while (1.0*rand()/RAND_MAX > y);
return x;
}
double random_theta() // Creating a random theta function
{
double x = acos((1.0*rand()/RAND_MAX - 0.5)*2.0); // Generates random No. between -1 and 1
return x; // acos() then translates into 0 -> pi
}
double random_phi() // Creating a random phi function, can call both later
{
double x = (1.0*rand()/RAND_MAX - 0.5)*8.0*atan(1.0); //Pi*random No. (-1 to 1) (i.e. phi range)
return x;
}
|