/* University College London Dept of Physics Course in C++ 3C59 
| all rights reserved 2000
|
| Main program to run the Compton experiment Monte-Carlo program
|
|
| Author: P.Clarke 
*/

#include <iostream>
#include <fstream>
#include <vector>

#include "Photon.h"
#include "Rod.h"
#include "PhotoMultiplier.h"


int main()
{

  // .................................................
  // Set some fixed numbers

  float distance = 50.0 ;        // distance of pmts from rod
  float incidentEnergy = 662.0 ; // Energy of source photons on rod
  float resolution     = 80.0  ; // Detector resolution in keV

  int numberOfEvents = 10000 ;     // Number of photons to generate

  // Ask for some input

  float electronMass ;
  std::cout << " Please input the electron mass in keV " << std::endl ;
  std::cin  >> electronMass ;

  float windowWidth ;
  std::cout << " Please input the width of the pmt window in cm " << std::endl ;
  std::cin  >> windowWidth ;

  // Create a Rod
  Rod myRod( incidentEnergy, electronMass ) ;

  // Create a vector of 18 photomultipliers at 10 degree difference from 0-180
  std::vector<PhotoMultiplier> detectors ;
  for( int ind=0; ind < 190; ind+=10 ) 
  {
    detectors.push_back( PhotoMultiplier( float(ind), distance, windowWidth, resolution ) );
  }

  // .............................................
  // Now run the experiment for many photons

  for( ind=0; ind < numberOfEvents ; ++ind ) 
  {
    // Create a new random photon
    Photon p = myRod.generatePhoton() ;

    // Offer it to each photomultiplier in turn
    for( int pmt=0; pmt < detectors.size() ; ++pmt )
    {
      detectors[pmt].tryit( p ) ;
    }
  }

  // ...........................................
  // Finally output the results by asking each pmt in turn for its average energy

  std::ofstream outfile( "mcCompton.dat" ) ;

  // Offer it to each photomultiplier in turn
  for( int pmt=0; pmt < detectors.size() ; ++pmt )
  {
    outfile << detectors[pmt].angle() << "  " 
            << detectors[pmt].averageEnergy() << "  "
            << resolution
            << std::endl ;

    std::cout << detectors[pmt].angle() << "  " 
            << detectors[pmt].averageEnergy() << "  "
            << resolution
            << std::endl ;
  }


  // Thats it 
  return 1 ;

}