/*
*          ANALYSIS OF BUBBLE CHAMBER PHOTOGRAPHS
*/

#include"lib351.h"

void main(){

    // output initialization
	ofstream dataOutStr;
    outFileMaker("m358p1c.out", &dataOutStr);
    
    // program variable definitions
	const double maxLTrack = 0.597;
	const double restMass = 1197;
	const double specificMomentum = 174;

	int eventNumber;
	double lTrack, th, s, lChord;
	double theta, curvMomentum, mMeVc, relativEn, lTrackMomentum;
	double xMass, relativEnTot;
	
	// write header
	dataOutStr<<"BUBBLE CHAMBER TRACKS RESULT"<< endl;

	// cycle through the runs
	bool goout = false;
	while (!goout){

		// data input
		cout <<"Please, supply:"<< endl;
		cout <<"\nidentifying number for the event, eventNumber:"<< endl;
		cin  >>eventNumber;
		cout <<"\nlength in cm of the Sigma track, lTrack:"<< endl;
		cin  >>lTrack;
		cout <<"\nangle in degrees between Sigma "<< endl;
		cout <<"& pi-minus track, th:"<<endl;
		cin  >>th;
		cout <<"\nsagitta in cm of the pi-minus track, s:"<< endl;
		cin  >>s;
		cout <<"\nchord length of the pi-minus track, lChord:"<< endl;
		cin  >>lChord;

		// write the data
        dataOutStr<<"eventNumber = "<<eventNumber<<endl;
        dataOutStr<<"length = "<<lTrack<<"; theta = "<<th<<"; sagitta = "<<s<<"; chord = "<<lChord<<endl<<endl;

		// calculate the result
   		theta = deg2rad(th);
		curvMomentum = pow(lChord, 2)/(8*s) + s/2;
		mMeVc = 6.86 * curvMomentum;
		relativEn = sqrt( pow(mMeVc, 2) + pow(140, 2));
		if((lTrack - maxLTrack) > 0){
			lTrack = maxLTrack;
		}
		lTrackMomentum = specificMomentum *
				pow((1 - lTrack/maxLTrack), (1./4));
		relativEnTot = sqrt( pow(lTrackMomentum, 2) + pow(restMass, 2));
		xMass = sqrt( pow((relativEnTot - relativEn), 2) -
	      	pow(lTrackMomentum, 2)- pow(mMeVc, 2) + 2 * lTrackMomentum
			* mMeVc * cos(theta));


		// write the answer
    	dataOutStr <<"mass of X particle = "<< xMass <<" MeV/(c*c)"<<endl<<endl;

		// another run?
		goout = query();
	}
    // close the output file
	fileClose("m358p1c.out", &dataOutStr);
};